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SUMMARY 


A numerical procedure is developed for the stress analysis of 
stiffened shells of arbitrary cross section. The procedure is 
based on the usual assumption of shell stress analysis—that the 
shape of the cross section is maintained by stiff bulkheads. The 
method is general enough to include the effects of cutouts, as 
well as the secondary effects due to restrained warping usually 
referred to as shear lag when the shell is loaded primarily in bend- 
ing, and bending stresses due to torsion when the shell is loaded 
primarily in torsion. 

The quantities determined directly are the spanwise distribu- 
tion of the axial displacements of the shell cross sections. The 
numerical equations defining these qualities are solved by sim- 
ple iteration. In order to obtain practical convergence of the 
process, it is necessary to refer displacements and loads to a set 
of axes in the cross section, which are defined as the principal 
axes for shear, and to resort to the occasional use of a special 
cycle in the iteration process, which in effect removes all displace- 
ments associated with planar distributions of axial displacement 
of the shell cross section. 

The details of the application of the procedure to a simple 
problem are shown, and the results obtained from a more com- 
plex problem are summarized. 


INTRODUCTION 


ee THE ORDINARY ENGINEERING THEORIES of 
stress analysis are used as approximate theories for 
the analysis of stiffened shell-type structures such as are 
used in aircraft structures, it is generally recognized 
that there are circumstances under which these theories 


can be considerably in error. Some of these circum- 
stances are more or less obvious. A concentrated load 
applied to a stringer cannot be expected to be immedi- 
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ately distributed throughout the rest of the stringers, 
or, if stringers are interrupted by a cutout, some con- 
centration of stress is expected in the neighborhood of 
the cutout. There are other circumstances that are not 
so readily apparent. Such situations are the shell in 
flexure when one end is built into a rigid support, 
when there is a spanwise variation of cross-sectional 
properties, or when the distribution of shear along the 
length of the beam is not uniform. Under these cir- 
cumstances the tendency of each cross section to warp 
out of its plane, which is necessary to establish the 
shear strains, is resisted by the rigid support or by the 
adjacent cross sections that want to warp differently 
and a change in stress distribution results. This 
problem has been recognized for some time and has be- 
come known as the shear lag problem. A similar 
situation exists in the case of the stiffened shell under 
torsion. If the shell is built in at one end or there is a 
variation in torque or cross-sectional properties along 
the length of the shell, the tendency of each cross 
section to warp is restrained and stresses are set up in 
the longitudinal directions. This problem also is well 
recognized, and the stresses induced have been called 
the bending stresses due to torsion. 

These two fundamentally similar problems have 
generally not been treated together in the literature. 
This, in part, has resulted from the fact that, in the 
development of those engineering solutions that are 
already available, a more simple physical picture of the 
problem has sufficed for the torsion problem than for 
the bending problem. Most of the solutions available 
have been limited to boxes of rectangular cross section 
in which the cross-sectional properties are symmetrical 
about both the horizontal and vertical axes. The ex- 
tension of these solutions to arbitrary cross sections has 
usually been too difficult to be of much practical use. 

The object of this paper is to establish a simple 
numerical method of stress analysis which requires but a 
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Fic. 1. Typical stiffened-shell structure of wing. 








Notation and coordinate system. 


Fic. 2. 


single point of view to handle the bending and torsion 
problem of stiffened shells. The analysis is based on 
the same simplifying assumptions as are found in the 
previously developed theories in this field but is capable 
of ready application to unsymmetrical cross sections 
and multicelled boxes. Because it was expected that 
convergence of the iterative process would not be rapid, 
the numerical process was set up so that at each step 
advantage could be taken of the process of continuous 
multiplication on an automatic calculating machine. 

The numerical problems actually solved are limited 
to the bending and torsion of single-celled boxes of 
rectangular cross section. It should be readily appar- 
ent that the method is general and can be applied to 
many other problems in stiffened shell analysis. 


THEORY 


Basic Assumptions 

The assumptions used in this method of analysis do 
not depart from those generally accepted for the engi- 
neering analysis of stiffened shell-type structures. The 
structure actually analyzed is a fictitious structure that 
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approximates the actual structure (see Fig. 1). The 
fundamental unit of the fictitious structure is a plane 
rectangular panel bounded on two opposite sides by 
rigid ribs and on the other two sides by stringers of 
finite extensional stiffness (see Fig. 2). These edge 
members are assumed rigid in lateral bending and are 
pin-connected at the corners. The panel itself consists 
of a sheet that has the fictitious property of being able 
to support only shear stress. The shear stress is con- 
stant within a given panel. 

The flexible stringers at the edge of the panel run 
parallel to the direction of the major stresses in the 
structure. For a panel loaded in tension the flexible 
stringers would be placed parallel to the tensile forces, 
Similarly, for a beam, these stringers would be placed 
in a direction parallel to the bending stresses. The rigid 
edges of the panels then lie at right angles to the major 
stresses and are in the direction of the stiffening nor- 
mally found in such structures. It is further assumed 
that these rigid members maintain the cross section 
and prevent any relative displacements of the points 
in the plane of a given cross section but offer no resist- 
ance to displacements out of the plane. 

With such a panel as has been described as the funda- 
mental unit, it is possible to build any type of stiffened 
structure as long as taper is excluded. The state of 
stress in such a structure can be completely defined 
when the following quantities are determined: 

(1) The displacements, at each end of each flexible 
stringer, in a direction parallel to the stringer. 

(2) The translation of the plane of each cross sec- 
tion defined by the rigid edge members. 

(3) The rotation of the plane of each cross section 
defined by the rigid edge members. 

The process of going from the actual structure to the 
fictitious structure is straightforward. The shape of 
the cross section is approximated by a polygon. AA side 
of the polygon may be the width of an entire panel of 
the actual structure or the sides may be subdivided 
into panels if they are large. The cross-sectional ma- 
terial resisting extension is then concentrated into the 
stringers that bound the edge of these panels. In 
general, the material in the half panel on either side of 
a stringer is included in the stringer area. This area 
may include real stringers, as well as sheet material. 
The thickness of the shear panels is taken equal to the 
thickness of the sheet. It would be more correct when 
choosing the fictitious stringer area to consider that the 
sheet material is more stiff than stringer material be- 
cause of the lateral restraint offered the sheet by the | 
bulkheads or ribs, but this was not done in the numeri- 
cal examples. 

Loads are applied to the structure only at the bulk- 
heads or ribs. All loadings can be considered either 
as axial loads applied directly to the stringers or as 
forces or twisting moments in the plane of the bulk- | 


head. 
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STRESS ANALYSIS OF 
SYMBOLS 


= width of jth panel 

i = subscript denoting bulkhead or bay, « = 0 at root and 
increases toward tip 

= subscript denoting stringer or panel location in cross 
section; proceeding clockwise when cross section is 

viewed from tip 


] = shift of origin in z direction when determining principal 
axis for shear 

m = shift of origin in y direction when determining principal 
axis for shear 

rj = length of perpendicular from center of rotation to 
panel 7 

;; = thickness of sheet in panel 7j 

uijj = axial displacement of stringer, positive when motion 
is away from root 

Au.” = in correction cycle axial translation of plane of ith 
bulkhead 

0; = displacement of ith bulkhead parallel to y axis 

Av; = change in displacement in 7th bay parallel to y axis 

w; = displacement of ith bulkhead parallel to z axis 

Aw; = change in displacement in ith bay parallel to z axis 

x = coordinate along length of box; increases from root 
to tip 

y = coordinate in cross section which corresponds to one of 
the principal axes for shear 

y”" = coordinate in cross section of bay which corresponds 
to one of the principal axes for bending 

Z = coordinate in cross section which corresponds to one 
of the principal axes for shear 

z" = coordinate in cross section of bay which corresponds 
to one of the principal axes for bending 

A;; = cross-sectional area of ij stringer 

Ap = total cross-sectional area of stringers in a bay 

E = Young’s modulus, arbitrarily chosen as 10 in illus- 
trative examples 

F;; = axial force applied to the stringers at the point 7j, 
positive when directed away from root 

G = shear modulus, arbitrarily chosen as 4 in the illustra- 
tive examples 

H; = total shear force in the #th bay parallel to the z axis 

I,” = principal moment of inertia of stringer areas in cross 
section about y” axis 

I“ = principal moment of inertia of stringer area in cross 
section about 2” axis 

eI = total shear stiffness of panels in a bay for displace- 


ments of bulkheads with fully restrained stringer 
displacements (see text for subscript notation) 
L; = length,of ith bay 


P;; = average force in stringer ij, positive when stringer is 
in tension . 

T; = total torque in ith bay about origin of principal shear 
axes 

V; = total shear force in the ith bay parallel to the y axis 

a = rotation of axis when determining principal axes for 
shear in cross section 

8,"; = rotation about y” axis of ith bulkhead for correction 
cycle 

8."; = rotation about s” axis of ith bulkhead for correction 
cycle 

vii = Shear strain in panel 7 
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Pitj =(6A) 
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2 P, =(@A); ; 
) | 
be 
2 NY 
Fic. 3. Forces on stringer. 
éij = average strain in stringer 7 
6 = rotation of ith bulkhead 
Ad; = change in rotation in the ith bay 
¢j = angle between positive y axis and rj; 
oij = average axial stress in stringer 7j 
tij = Shear stress in panel ij 
2; = summation around periphery of cross section 


STRAIN RELATIONS 


The following relationships will be assumed to exist 
between the above-mentioned displacements and the 
strains and stresses. (See “Symbols” and Fig. 2 for 
the notation and coordinate system used.) The shear 
strain in the 77 panel is assumed constant and given by 
the relation 


¥9 = (lag 4a + M4154) - 
r; 


1 A Av, 
(ui,g + Ui 41,9)] 2%, + on ae > 


and the shear stress by the relation 
(2) 


Because of the constant shear stresses in each panel, the 
normal stresses in each stringer will vary linearly, and 
the average strain in the stringer is given by the ex- 


Ti, , = Gy: r] 


pression 
€,9 = (Ui 41,3 — Me, 9)/Lt (3) 
and the average stress by 
o,,= EG; (4) 


EQUILIBRIUM OF STRINGERS 


Consider that portion of the structure immediately adjacent to the point 7j and isolate from it the stringer lying 
half a bay length to either side of the point ij. The forces acting on this portion of the stringer are shown in Fig. 3. 
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That the average force exists in the stringer at the center of the bay follows from the fact that the stresses in each | 
stringer vary linearly. 
The static equilibrium of this portion of the structure yields the following equations: 


- - L; —1 
(Poy — Peas) — (Oe - 1.5 1 — (Os 1, 4] alla [(rt)i,3 —1 — (rt)s, Az “+ Fy = 0 (5) 
Substitution of the following relations in Eq. (5) yields a relationship among w,,, the displacements at the ends 4 
of the adjacent stringer section, and the displacements of the ij cross section relative to the two adjacent cross 
sections. 
P,, = (@A)y = (Ui 41,5 — 4, )(AE/L), 3 
Ni, ji 441, paca) (tes 4 + Us +1, 9) rj Av; = Aw; 
(rt) = (Gt); — —-+ Ad, — — sin ¢; + — cos¢ 
[' 2b, 2b, ae * as * 7 4 
This relationship is given by the expression 





sf GtL GtL 
Ui, 4 = Fs" pee Mga —s + (tes, 94-1 M041, 5 +1) esx; 
| (4 z) (‘ i) (=) 4 | (“ “) (2) (=) | r A 
- ~ — u 7 — — 
idinia L/i-1, tb J; ~1, 1b J -15-1 one L/i3 4b Js, j 4b J i,4-1 ch 
‘ (a2) ied ‘. Cs) - C3) B to 
“ui 1, j 1 Uj, j 1 th ee Ui, j-1 u it, paokl 1b - i-1l 9 poet = 
Gtr Gtr Gtr Gt sin ¢ Gt sin ¢ 
— + Ad; — — Av;-1 — 
Z im}, f 1 2 i, j 2 i,j-l i—i, j 2 i—1l, j—1 
a» (“ = *) - (= = “) | oo es cos “) (2 = j | 4 
2 if 2 t, 1 ye “ iat, 93 
((A a) - (“8 cos *) | 4 F, (6 
i gat 


where in| 
rs = |(“ “) Ps (“ =) 4 (=) in (=) 4 (Sy) 4 (=) | 
} L i—1, j L ae tb i—l, j—1 4b i-—1, j 4b i, j—1 4b ij 
The above relationship assumes no equality of any of the structural elements surrounding the point 77. For | : 
the structure with constant cross section and bays of equal length, the expression for u,; reduces to g tor 
) the 


"| GtL Gtr Gtr 
(t6¢—1, ga + 26s, ge H M141, s-1) tb + (40,1 + Ad,) —) — (— « tink ees 
: = J - j-1 
Gt sin ¢ Gt sin ¢ Gt cose i ean 
(ee ae eee oe co 
2 j 2 , 5 9 - 


P L Su ome “ (S) + (u re (4 F) (=) (=) | | will 
=< i j+ oui, j41 itl, f+ i—1,j ¢i +1, + ats af iowa re 
t, J vs 1, §+1 1, f+ 1b 1 Se a 4b /, 4b), 1 ¥ ; os 





Where some of the structural elements are lacking because of cutouts, the stiffness associated with these elements I 
go to zero and may eliminate the corresponding displacements from the above expression. give 
the 
STATIC EQUILIBRIUM OF THE Bay : 
- 
° ° oe ° § Joe 
In the expression for u,;, the quantities A@,, Av;,, and Aw, are required. These quantities can be obtained by | 
relating the shear stresses in each panel to the total torque and shear forces acting in that bay. 
If T,, V;, and H, are the torque and shears in the 7th bay, then A 
’ Je 
T; = DL (rt):, jor, Vix —Di(rt);, jo; sing, H, = DY(th);, 6; cos ¢; 
j J j 
Substitution of the relationship between the shear stress and the displacements of the structure yields the following Je 
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Gtbr? Gtbr sin ¢ Gtbr cos ¢ 
DT) ~ LAE), + (9), - 
ee eee L o/s ‘hon EL ses 
j j j 
- > | Gtr Gir 
T+ (ui4s, 4 + Uj, a (“) = (=) | 
</i3 “ /t, 4-1 
Jj 
Gtbr sin ‘) } (2 sin? *) Y (F sin ¢ cos ‘) 
6; = - — Av, — Aw ——_——— = 
‘ »( iL i, J ; . pa i, Jj > ; . 8 
J J Jj 
Y : | (Gt sin ¢ Gt sin ¢ 
“> > (Mi4s, j tu, | (FE *) = ( 7 *) ] 
7 & i, J “ t, j-—1 
J 
Gtbr cos *) (= sin ¢ cos ‘) (“ cos? *) 
6 ————*) — wy, ) — ‘ Aw > —— = 
wd ( L i, j L ij * " : L i, j 
7 J 
Gt cos ¢ Gt cos ¢ 
no Door nf), (224) J 
- - i, j - 1, j-1 


A consideration of these three equations of statics indicates that they can be considerably simplified by the proper 
choice of the origin and orientation of the axes in each bay. By a proper choice of axes (see next section) the inter- 
action terms among the rotations and displacements can be made to disappear so that the same relationships reduce 


to 
l Y ; Gt 
I *te + (iy + Ui+1, A(Z 
/ 6,6 r 2 




















dé; 


Bad | 
tj 2 i, j-—1 

re Gt sin ¢ Gt sin ¢ l 
Av; = (a it " ees Does + Ui41, al ( 2 f) , — ( r *) a} (8) 

j 

1 3 : Gt cos ¢ Gt cos ¢ 

Aw; = Len + (us, 3 + Uj +1, l(‘ 2 *) = (¢ <- *) } 
J 00 ; - i, J - t, j—1 ) 


Gtbr? } Gtb sin? ¢ , } Gtb cos* ¢ 
Joo = >; Joo = (= “ae ) > Juw = ( : = ) 
- > ( L ) 3 a  & i, j ‘ 7 L i,j 
j j j 


The new set of axes will be referred to as the principal shear axes. All cross-section displacements, section 
torque and shears, and stiffnesses will be referred to these axes. A simple physical significance is associated with 
these axes. If the stringers in a given bay are rigid and held fixed at one end, a torque at the other bulkhead 
will produce only rotation about the principal shear axes, while a shear load in line with one of these axes produces 
only a displacement in the direction of the axis. 








in which 





DETERMINATION OF THE PRINCIPAL SHEAR AXES i } (2 sin yg’ cos e) 
°.w = a = — 
i 


In order to determine the principal axes for shear in a j L 


given bay, select an arbitrary set of axes and compute 


ee in which the primes refer to the original choice of axes. 
the quantities 


Substitution of the relationships (see Fig. 4) 


f- 2 atbr’ si : , . 
hij ) (F ‘) > Jy,’ = ) (+ oi . ¢ =e¢gta; r =r+Ilsing’ + moos gy’ 
ls , , j 
j ad Jj 


in the above stiffness formulas and setting 


, , 7 } Gtb sin? g’ 
Jou’ = (Gtbr’ cos ¢’);3 dig = ( ao 2) Jo» = Jo = Jo = 0 
. - 4 


j ‘ for the principal set of axes, yields the following rela- 


; (2 cos? “) tionships: 
bos: tan 2a = 2Jeu"/(Juw’ — Ing’) 
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/ ; / , , 
l SF tg J 0.0 aes J 6, Jou J 6,0 Fes _- Jo» Jom 
= . m = 
, , 9 ’ , 9 
ay Fos = Se © Fas Jou ae Jun” 
Joo = Joe’ — lJo. — mJ o. 
, , , , 
Je eee es — Poe 
Pe + (cos 2a) — 
» ») 
Faas (sin 2a) 
and 
, , , 
Jae + Jan Jee — Jum 
J eo _ (cos 2a) + 
Jy’ (sin 2a) 


CORRECTION CYCLE 


Although the sets of Eqs. (6) and (S) are sufficient 
to define the displacements of the structure, it was found 
that in trying to iterate these equations the process of 
convergence was, in general, impractically slow. It 
was also apparent that the difficulty was associated 
with displacements Av, and Aw; To speed up the 
process a correction cycle was developed which proved 
extremely effective. 

The basis of the correction cycle is the fact that in 
any bay the stringer forces must be in equilibrium with 
the thrust and moments in that bay. The correction 
cycle consists of adding to the strains in any bay a set of 
strains which is planar-distributed in the cross section 
such that the condition of thrust and moment equilib- 
rium is established. The displacements Av; and Aw, 
are then adjusted to restore shear equilibrium of the 
For the shell loaded in pure torsion, there is no 


bays. 
This is also 


resultant force or resulting moments. 
true if the solution sought is a correction to be added to 
a solution that already satisfies static conditions but 
not compatibility of strains. The conditions of equilib- 
rium on the cross sections for such solutions are: 
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Ps = 0; UP 7 = 9; UP as 2,5 = 0 (1) 
J J J 
The set of correction strains times bay length can be 


represented by the quantity 
Au; + Bai; + By Bj 


in which Aw, is the relative axial translation of adjacent 
bulkheads and 8., and £,; are the relative rotations 
about the axis zs and y of the planes of the bulkheads, 
Substitution of this type of displacement in Eq. (11), 
along with the set of displacements to be corrected, 
vields 


oi 


AE 
Do (us 43159 — May + Aus + Bad) + re) ) = 0 


AFy 
Dd (ui 21,9 —~ 4,3 + Au; + Bi; + Bss)( I *) = 0) 
7 4 


Do (uigy — My + Aus + Bays + ane)(“) = 0 
j L i, j 
Consideration of these equations indicates that, if the 
coordinates y and z are measured from the principal 
axes for bending of the cross section, a considerable 
simplification can be made for the determination of 
Any, Be, and B,;. Making use of the principal axes for 
bending and dropping £ and L yields 


i 
] . 

Au; = - \ > (a2, 9 u;, (A)s, 3 
447 3 

i I f ” \ 9 

B."; = “< Do (uigs, y — Ue (Ay")s, 5 (12) 
yij 

3.” s a a 4-" 

Bye = ~F7 Leia om Me JAS as | 


Se 


in which A7 is the total stringer cross-sectional area 
and J,-,; and J,; are the two principal moments of in- 
ertia of the 7th cross section. Because the above rela- 
tionships determine only the difference in displacement 
between bulkheads, there still exists the problem of 
how to add these to the original set of displacements. 
This is done by starting in the bay at the fixed end. 
When the set of corrected values for the first bay 1s 
obtained, they are added to the displacements at the 
first bulkhead beyond the root. With the corrected 
set of displacements for the first bulkhead, a set of 
corrections is determined for the second bay. The first 
bulkhead is considered fixed, and all the corrections are 
added to the bulkhead. This procedure 1s 
continued to the tip. After the corrections have been 
applied to the stringer displacements, the first equa- 
tions iterated are these for A@;, Av;, and Aw,. There 
should be no change in Ad. 

The following physical interpretation describes the 
action of the correction cycle. The planar changes of 
the stringer displacements removes the unbalanced 
moments and thrusts that have developed on the cross 
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sections. In doing so, they have given rise to changes 
in the shear strains in the individual panels. It can be 
shown that these changes in shear strain correspond to 
unbalanced shear forces on the cross section but no un- 
balanced torque. The following iteration of A@,, Av, 
Aw, will therefore indicate no change in AQ, as well as 
the necessary change in Av, and Aw, to remove the un- 
balanced shears. It can also be shown that there is no 
change in the shear strain in any panel due to the appli- 
cation of the correction cycle. 


NUMERICAL EXAMPLE 


In order to demonstrate the application of the numeri- 
cal method of analysis, consider the simple box beam 
shown in Fig. 5. It is a simple rectangular box sym- 
metrical about a horizontal axis. The corner flanges 
are not of equal area, so that the box is unsymmetrical 
about a vertical axis. The analysis is set up for two 
types of loading: a concentrated torque at the tip and 
a concentrated vertical load applied at the elementary 
“shear center.”’ According to elementary theory of 
bending, the vertical load should produce no twist of the 
box. 

Because of the symmetry of the box about a hori- 
zontal axis the displacement of corresponding stringers 
in top and bottom of the box are equal and opposite for 
Therefore only one-half of 
the box need be considered. Although the stringer 
material is unsymmetrical about a vertical axis, the 
shear webs are symmetrical about both a horizontal and 
vertical axis, and the principal shear axes coincide 
with the centerlines of the box. The bay shear stiff- 
nesses are calculated to be 


the loadings considered. 


Joe = >>(Gtbr?/L) = 70.0 


J 


Joo = Do(Gtb sin? g/L) = 0.2667 
j 


The horizontal shear stiffness of the bay will not be 
required for the loadings considered. The stringer 
stiffnesses are 


(< “) (“ =| i 
= = ().6667; = 0.1667 
ae ee L /:,2,3 


(AE/L),,4 = 0.3333 


The shear panel stiffness for axial displacements are 


GtL GtL 
2 = ().2000; = (0.1000 
4b i 4b i, 2 


(GtL/4b) ;, 0,4 = 0.3000 


The shear panel stiffnesses for rotation of the bulkheads 


are 
Gtr : Gtr 
- = 0.5000; (—> = 3.000 
© 741,23 ~ /4,0,4 


The shear panel stiffnesses for vertical translation of the 
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Fic. 5. Details of idealized structure, 


bulkheads are 


Gt sin ¢ Gt sin 
( : *) = —0.2000; ( ‘) =0 
2 i, 0 2 i, 1,2, 3 


(Gt sin ¢/2);,4 = 0.2000 


The correction cycle consists of only a rotation about 
the horizontal axis. This follows from the fact that 
symmetry automatically satisfies equilibrium of axial 
thrust and moment about the vertical axis. The prin- 
cipal moment of inertia about the horizontal axis for 
half the box is 


I, = (1/2))5(Az");, 5 = 100.0 

7 
Substitution of the computed stiffnesses in the general 
equation of equilibrium for the stringers yields the 
following set of equations: 


Us,1 = [(ui—1,2 + 2us,2 + u145, 2)(0.6818) — 
(4— 1,1 + Ui+1, 1) (0.04545) — (A0@;-1 of Ad;) x 
(0.8523) —(Av;_1 + Av,)(0.06818)] 
Ui,2 = [(uj—1, 3 + 2u;3 + Ui41,3)(0.1071) — 
(uy-1,2 + Mi41, 2)(0.1429) + 
(ul; iS + 2u4,1 + Uj+1, 1)(0.2143)] 
Us,3 = [(ui—1,4 + 2m 4 + ui 41, 4)(0.2143) — 
(Ui~1,3 + Ui+1, 3) (0.1429) a 
(up—1,2 + 2,2 + W441, 2)(0.1071)] 
“14 = [(u,- 1,3 + 2:3 + U;41,3)(0.08824) — 


(2;- ~ a Uz 41, 4)(0.2059) —f- ( A@;- it Aé,) X 
(1.1029) — (Av;_1 + Av,)(0.08824)] 


The equations for u;; and u;,4 include the effect of 
assuming that the flanges on the lower side of the box 
have equal and opposite displacements. The equations 
apply at points in the interior. Special equations could 
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be written for the free ends of the stringers at the tip 
and root. The above equations, however, will also 
suffice at the tip and root if it is imagined that there 
exists beyond the ends a fictitious bulkhead and bay 
that have the same displacements as the first bulkhead 
and bay inside the ends. 


The equilibrium equations for the bays are for bend- 
ing, 
Ad; = —38.3 + (0.07143) (—u, 1 @+1,1 + 
U4,4 + Ui +1, 4) 
Av, = —3750 — (1.500) (ui. + ue+i,1 + 


Ui,4 + U; +1, 4) 


and, for torsion, 


Ad, = 14.38 + (0.07148)(—u41 — upti.1 + 
Ui4 + Uy +1, 4) 
Av, = —(1.500) (ui. + ug4i.a + Ue + i 41,4) 


The details of the iteration for the box loaded in bending 
when only the flanges are attached at the root is given 
in Fig. 6. The original assumption for stringer dis- 
placement is given in vertical lettering in the upper 
left-hand corner of the box representing each point 
and corresponds to the engineering theory of bending 
displacements. The original assumptions for the bay 
displacements were obtained from the original choice 
for the stringer displacements by application of the bay 
equilibrium equations. The values just below the 
original assumption were obtained by iterating the sets 
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Details of the iteration for bending of idealized structure with partially free root. 
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*Correction cycle 


(Final values are boxed 


of equations, using only the original set of displace- 
ments. The following number, which is slanted, is the 
difference between the original displacement and the 
first iterated displacement. Succeeding values are the 
successive approximations for the difference between 
the true solution and the original assumption. 


original set, except that any constant terms are elimi- 
nated and are replaced by the difference between the 
original assumption and the first iterated value. For 
instance, the equation for the difference solution for 
Ui, 3 is 

Ui.3 = —361 + [(ui—1,4 + 2us,4 + i +1,4) (0.2148) - 


1,3 + uj t 1, 3)(0.1429) + (ui —1, 2 + 2uj, 2 + 
U; +1, 2)(0.1071 


(uj 


and for the vertical bay displacement 


Av; = —(1.500) (uy. + upti,1 $a + Ms t14) 


The principal reason for iterating the difference 


solution is that the difference solutions correspond to 


self-equilibrating state of stress which makes it simplet 
to apply the correction cycle. It also has the advat- 
tage that smaller numbers are required in the succeed 
ing iterations. 

Iteration was started at the point m2 and continued 
around the root bulkhead before moving to bulkheaé 


one. After the stringer displacements were iterated, 
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Fic. 7. Distribution of stringer stresses (pounds per square 


inch) for bending of idealized structure with partially free root 
for four and infinite number of bays in length of 120 in. 


improved values until the sixth cycle, at which time a 
correction cycle was applied. 

If the root had been fixed, the correction cycle would 
have started at the first bulkhead assuming the root 
fixed. However, because some of the stringers at the 
root are free and improved displacement is available 
in bay zero and at bulkhead one, these displacements 
at the root are iterated normally before the root is con- 
sidered fixed and before the corrections are applied to 
the outboard bulkheads. The first set of corrections 
is obtained by applying Eq. (12) to the first bay: 

— J 
8,"y = ~ [(563 — 0) 5 K 2 + (—273 — 3979) 5 X 
100.0 
0.5 + (—124 — 4102) 5 X 0.5 + (1073 — 0)5 X 1] 
5 = 510 


This correction term is applied to the displacements at 
the first bulkhead to yield the corrected displacements. 
In a similar fashion the entire set of displacements are 
adjusted. The adjusted displacements are given by the 
first set of starred numbers. The adjusted bay dis- 
placements were obtained from the adjusted stringer 
displacements. 

The final accepted values of the differences solution 
are given as the last slant figure in each box. The final 
accepted displacements, which are the sum of the dif- 
lerence displacements and the originally assumed dis- 
' placements, are enclosed in the small boxes at the end 








& 


of each list of iterated values. 

Stringer stresses corresponding to the sets of dis- 
placements just obtained for this loading condition are 
For purposes of comparison, solutions 
obtained for this same structure and loading but con- 
sidering an infinite number of bays are also included. 


rhese solutions are given in Appendix A. Fig. 8 shows 
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Fic. 8. Distribution of stringer stresses (pounds per square 
inch) for torsion of idealized structure with fixed root for four 
and infinite number of bays in length of 120 in. 
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plots of stringer stresses obtained by this same method 
when the tip is loaded by a concentrated torque. 
There is also included in the figure the solutions ob- 
tained for an infinite number of bays. 

The previously described numerical example was con- 
siderably simplified in order to demonstrate the method 
An example, which is a more critical test of 
the application of the method, is the stress analysis of 


involved. 


the structure shown in Fig. 9. No details of the analy 


sis are shown, but a summary of results is given in Figs. 
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10 and 11. In Fig. 10 a spanwise plot of stringer 
stresses is given for this box loaded in torsion at the tip 
and fixed at the root. Stresses are given when two and 
four bays were assumed to exist along the length of the 
box. In Fig. 11 the final stresses are shown for the 
same box loaded by a concentrated force at the tip 
when only the four corner flanges are held fixed at the 
root. 


CONCLUDING REMARKS 


The numerical procedure that has been described 
briefly has been applied to other structures and loadings 
than the ones considered here. These are the plane 
stiffened panel: loaded in bending or tension by con- 
centrated loads in its plane and the plane panel, with a 
large cutout, loaded in shear or in tension. It has also 
been used to compute secondary stresses induced by 


nonuniform temperature distributions. Satisfactory 
convergence of the iteration process has been found in 
all these cases. 

It should be apparent that the application of the 
method to multicelled boxes or to open-celled boxes 
can be made without difficulty. 
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Appendix A 


DIFFERENTIAL EQUATIONS FOR THE UNIFORM Box BEAM WITH AN INFINITE NUMBER OF Bays 


In order to obtain the solution for the idealized structure considered in the numerical example when the bay 


length is zero, the differential equations for the problem were developed and solved. 


The differential equations 


are readily obtained from the equations already developed for finite bay lengths. 


If Eq (7) for the structure of constant cross section is rearranged and divided by the bay length, the following 


ee | (2) ‘ (S) |+ 
eS ea Se 


expression is obtained: 
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If L is allowed to approach zero, the following differential equation results: 


pn | (2) +(2)] ~ (“) +(%) + [Gin — Grr), 1S 
(AE); 73 Die b/, uy; b aig b mee r); (Gir);—1) 7 









































E dv dw dF 
[(Gt sin ¢); — (Gt sin ¢); -|J— + [(Gt cos ¢); — (Gt cos ¢);-,] ~——! = @ (A-1) 
dx dx dx 
Similarly, Eqs. 8 can be reduced to the following form: 
do/dx = [1/>°(Gtbr?),] {T + You,[(Gtr), — (Gtr); -1]} 
j j 
‘yg dv/dx = [1/>5(Gtb sin? ¢),] { V — Yu, [(Gt sin ¢); — (Gt sin ¢);-,]} (A-2) 
j j 
dw/dx = [1/>°(Gtb cos? ¢),] {H + Sou,[(Gt cos ¢); — (Gt cos ¢);-1]} 
j J 
The expressions for the stringer and shear stresses become 
o,; = Edu,/dx 
r square ) 10 F F 
root f | (@j+1 — 4 t - av aw : A- 
7 1 = d at Ss vj — — sin yg; — + cos 3] (A-3) 
b; dx dx dx 
factory An interesting result from Eqs. (A-2) and the expression for the shear stress is that, at a completely restrained 
und in | section where wu is zero, the shear stresses are completely defined by the external forces and the shear stiffness of the 
cross section. 
of the 
| boxes FF THE TORSION PROBLEM WITH FIXED ROOT 
By substituting the numerical values of the cross-sectional properties and a tip torque of 1,000 in. Ibs. into the 
differential Eqs. (A-1) and (A-2), the following simultaneous differential equations are obtained : 
' (d2u,/dx?) — 0.0021429%, + 0.001333u. + 0.0008095u, = 0.11905 
ciation | (d*u2/dx?) — 0.008000u2 + 0.005333, + 0.0026667u;3 = 0 
f Civil (d?u3/dx*) — 0.008000u; + 0.0026667u2 + 0.005333u = 0 
sine (d?u4/dx?) — 0.0042859u, + 0.0026667u3 + 0.00161905u, = 0.23810 
e been The solution of these equations for the conditions of fixed root and the assumption of an extremely long 
uks to span length is 
andard ; se ; 
ttitude U, = 1.365e—°-1 — 1 552e—0- SH + 39. 243e—0.080% — 39.057 
us = —8.064e—%. 111% +. 8 230e—9.85% + 12. 849e—9.0589 — 13.015 
Uz = 9.822e—9. 1119 4. 5 OSle—-0859 — 53. 92Ze—0.058% +. 39.069 
4 = —3.613e-9-11% — 3, 52G6e-%.999% — 57. 967e-9-95%% 4+ 65.110 


The corresponding stringer stresses are 


oy = —1.527e—%-119% 4 1 333e—0.0850% 4 21 14 1¢e—0.053% 
= 9. 020e-9-1119 — 7. O7D0e—9-9859 — §. 922e—0.05302 
— 10. 986e~%- 1119 — 4,322¢—0.0859% + 29 048e—0.058%2 
oy = 4.041le—%- 1182 +. 3 029e—9. 0859 + 3]. 227e—0. 05302 


2 
| 


q 
nw t 
ll 


the bay 
uations 


The numerical results obtained from the equations for the stringer forces are plotted in Fig. 8. 





Llowing 
THE BENDING OF THE BOX BEAM WITH PARTIALLY FREE ROOT 

i For the bending problem it will be assumed that there is at the tip of the box beam, a force of 1,000 Ibs. applied 

f ° . ° ° ° 

t at the shear center. The differential equations for this case are the same as for the torsion problem, except the 

right-hand sides now become 
. + Y —2.8188 
0 

| + 0 
1 —4,3622 


The following results for the stringer displacements were obtained when it is assumed the box is long and only 
the corner stringers are attached: 
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uy, = 300x — 1.25x? + 69.0e—9-1119 — 773. 7e—9-915 — 390. 9e—9-998 + 1,095.6 

us = 300x — 1.25x? — 407.8e—9- 1! + 4,104. 0e—?-9359 — 128. 0e—?-9% + 693.8 
uz = 300x — 1.25x? + 496. 7e—9- 1! + 2,508. 5e—9-9859% 4+ 537. 2e—9-05897 + 897.7 
us = 300x — 1.25x? — 182. 7e—9-1119 — 1,758. 2e—0. 08592 4. 577. 5e—0.0539% 4. 1 363.5 


The stringer stresses corresponding to these displacements are 


o, = 3,000 — 25x — 77. 2e~%-!% 4 664. 6e-%-" + 210. Ge-0.098 
o2 = 3,000 — 5x + 456.0278." — 3,525.02-0.m¥* + 69.06 0,052 
a3 = 3,000 — 25x — 555.6e--1% — 2,514. if — 299. 4e—0.05802 
os = 3,000 — 25x + 204.3e-%11% 4 1,510. 2e-.985% — 311, 1e—0.0589% 


These stringer stresses are plotted in Fig. 7 
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Two-Dimensional Jet Mixing of a 
Compressible Fluid’ 


S. I. PAIt 


Cornell Aeronautical Laboratory, Inc. 


SUMMARY 


The mixing and divergence of a supersonic jet exhausting into 
a supersonic stream are investigated theoretically. 

In the first part of this paper, the flow is assumed to be laminar. 
When the velocity and temperature in the jet are different 
slightly from those of the surrounding stream, by the method of 
small perturbations and under ordinary boundary layer assump- 
tions, the equation of motion of two-dimensional flow will be 
reduced to a form of the well-known equation of heat conduction, 
whose solution is known for any given boundary conditions. It 
has also been shown that the exact solution of the two-dimen- 
sional jet mixing of viscous compressible fluids can be obtained 
by successive approximations starting with the solution of small 
perturbations. 

Velocity and temperature distributions for two cases—one is 
the mixing of two uniform flows and the other is the mixing of a 
jet of compressible fluid from a two-dimensional nozzle with full 
expansion exhausting into a supersonic stream—have been cal- 
culated. The properties of the jet mixing depend mainly on the 
momentum of the jet regardless of whether the change of momen- 
tum is due to the change of velocity or the change of tempera 
ture—i.e., the change of density. Compressibility has a con- 
siderable effect on the properties of the jet. 

In the second part, the cases of turbulent flow are investigated. 
By means of Reichardt’s theory of free turbulence, the turbulent 
shearing stress may be expressed as 


tT = pe(Ou/Oy) 
and 
€ = k(tmar. — Umin-)b 
It has been shown in this paper that 
«= €9(Xx, L)* 


where ¢9 is a constant that can be determined experimentally. 
The value of n lies between 0 and 1. The exact value of nm de 
pends on the condition of mixing. 

When the expression of turbulent shearing stress given above 
is used instead of the viscous stress in the equation of motion, 
by suitable transformation of variables, it has been shown that 
the equation of two-dimensional turbulent jet mixing is identical 
to that of the laminar case. Hence, the solution of the first part 
of this paper can be applied to the turbulent case, provided that 
the characteristic constants ¢9 and m have been properly chosen. 


List OF SYMBOLS 


A = constant 
b = width of the mixing zone 


Presented at the Fluid Mechanics Session, Seventeenth Annual 
Meeting, I.A.S., New York, January 24-27, 1949. 
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theoretical interest and practical importance. 
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B = constant 
Cp = specific heat at constant pressure 
Cy = specific heat at constant volume 
f = p*p* 
fi = u,*yn*p* 
F = integral defined by Eq. (21) 
g = gravitational acceleration 
G = integral defined by Eq. (23) 
k = constant 
L = a characteristics length 
m = constant 
n = constant 
a = pressure 
R = gas constant 
T = absolute temperature 
To = absolute temperature of the surrounding stream of 
velocity uo 
Tr = Absolute temperature of the flow at velocity u, 
u = x-component of the velocity 
uo = u of the surrounding stream 
ah = u — Uo 
ur = largest u in the jet 
1410 = ul — Uo 
a* = u/Uo 
u,* = u,/Uo 
Umar. = Maximum 4 at a given station of x 
Umins = Minimum uz at a given station of x 
v = y-component of the velocity 
v = y-component of the velocity 
£ = coordinate taken along the axis of the jet 
Bg = coordinate defined by Eq. (35) 
y = coordinate taken along axis perpendicular to the 
x-axis 
vo = value of yat 7 = 0 
pated 
a = V po/pouo 
€ = coefficient of eddy kinematic viscosity 
€0 = constant having the same dimension as € 
n = variable defined by Eq. (15) 
m = value of n atx = Oand y = 1 
= coefficient of heat conduction 
u = coefficient of viscosity 
Ho = coefficient of viscosity at 7 = 7% 
u* = u/uo 
p = density 
po = density of the fluid at T = 7) 
pr = density of the fluid at T = 7, 
p* = p/po 
T = turbulent shearing stress 
@(y) = (2/V w) So'e-® dt = Gauss’s error integral 
y¥ = stream function 
v1 = value of y at x = Oand y = 1. 


INTRODUCTION 


4 I ‘iE PROBLEM OF THE FLOW of a jet of a compressible 


fluid exhausting into a uniform stream is of great 
In this 
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paper both the laminar and the turbulent flow of two- 
dimensional jet mixing of a compressible fluid are in- 
vestigated theoretically. 

In the case of high-speed flow from a jet engine, the 
density of the fluid is low and its temperature is usually 
high so that the kinematic viscosity is large; hence, the 
resulting Reynolds Number will be relatively small in 
spite of the high speed. At relatively small Reynolds 
Number, it is expected that a considerable portion of the 
jet would be laminar. Hence, the solution assuming 
laminar flow would be of some practical interest. As 
will be pointed out in this paper, the results for laminar 
flow are applicable to turbulent flow provided that the 
empirical quantity of eddy kinematic viscosity has been 
experimentally evaluated. This is another reason for 
investigating the laminar jet. 

Turbulent flow is usually important in most prac- 
tical cases, particularly far downstream. However, 
our knowledge concerning the laws of turbulent flow of 
a compressible fluid at high speeds is still in a lamentable 
We have to use the so-called semiempirical 
But as far as the 


state. 
theory to analyze the turbulent flow. 
mean velocity distribution is concerned, such theories 
usually give fairly good results. In this paper, Reich- 
ardt’s theory of free turbulence has been extended to the 
case of a compressible fluid. As usual for the semi- 
empirical theory, it is necessary in the final results to 
evaluate experimentally some empirical characteristic 
constants, such as the coefficient of eddy kinematic 


viscosity. 
THE LAMINAR JET MIXING 


Fundamental Equations 


The equation of motion of two-dimensional jet mix- 
ing of a viscous compressible fluid is 


Ou _ Ou o ot) 
re Oe” a (. ay 

This equation is obtained under the ordinary assump- 
tions of boundary layer and the pressure is assumed 
to be constant. The axis of x is taken along the axis 
of the jet and the axis of y perpendicular thereto. The 
quantities w and v are the x- and y-components of the 
The density as designated by 


(1) 


velocity, respectively. 
p and pis the coefficient of viscosity of the fluid. 
The equation of continuity is 


[O(pu) /Ox] + [O(pv)/Ov] = O (2) 
The equation of energy is 
o Cp T) + oO (CT >(r or) * (**) 
Uu : > ) v L)= / " 
a . , Ov ‘ Oyv\ Ov ‘ Oy 
(3) 


where Cp is the specific heat at constant pressure, 7° is 
the absolute temperature, and d is the coefficient of heat 
conduction. 
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Flow of a two-dimensional jet in a uniform stream of 
velocity Up. 


Fic. 1. 


It was first shown by Crocco® ? that, if the Prandtl’s 
number Cpu/X is assumed to be unity, the following re- 
lation satisfies Eqs. (1) and (3) 

T =A + Bu — (u?/2Cp) (4) 


where A and B are constants determined by the bound- 
ary conditions. 
The pressure is assumed to be constant; hence, 


p* alt. ‘po - To/T = 1/T* (5) 
The coefficient of viscosity may be expressed as 
u* = w/o = (T/To)™ = T*” (6) 


where m is taken as a value less than or equal to one.” * 


Solution by the Method of Small Perturbations 


In some practical cases, both the velocity and the 
temperature in the jet are only different slightly from 
those of the surrounding stream (Fig. 1). It is ad- 
visable to find some approximate solution for this spe- 
cial case first. 

Write 


u= uUu+ M1, v= (7) 


where #1, << up and u ~ 2. 
Substituting Eq. (7) into Eq. (1) and neglecting the 
higher order terms, we have 


Om,/Ox = a*(0°u,/Oy") (8) 
where a? = puo/poto. 

Eq. (8) is the well-known equation of heat conduc- 
tion whose solution is known for any given boundary 


conditions. 4 


Mixing of Flow from Two-Dimensional Nozzle Under Full 
Expansion 


from a two-dimensional nozzle 


i.e., the pressure of the flow at 


Consider the flow 
under full expansion 
the exit of the nozzle is exactly equal to that of the 


surrounding stream. Then, the velocity distribution 
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at the mouth of the jet, Fig. 2, may be assumed as fol- 


lows: 


5 = & | 


u = Uw, —l< y< ] 
_1\ 


1< yandy< 
The solution of Eq. (8) with the boundary conditions 


(9) 


u, = 0, 


of Eqs. (9) is 


_ Mo E (; = *) 14 ( +2) | (10) 
2 QaVx 2aVx 

» "9 be 

~ | e *dé 

Vare? 


equals Gauss’s error integral. The result of Eq. (10) is 
The divergence of the jet is propor- 


where 

o(y) = 
plotted in Fig. 3. 
tional to a = V po/ potto. 


Mixing of Two Uniform Flows 


The velocity distribution at the beginning of the 


mixing zone (Fig. 4) is 


x =0 
t, = 0, y >0? (11) 
uy = Uy, y< of 


The solution of Eq. (8) with the boundary condi- 


tions of Eqs. (11) is 
“ 
WM = (uw/2)[1 — o(y/2aV x)] (12) 


The results of Eq. (12) are plotted in Fig. 5. 


Exact Solution 


To solve the problem more rigorously, one has to re- 
sort to Eqs. (1) and (2). By introducing the stream 
function y, which is defined by 

OY OY ' 
= pas = —p*y (13) 
oy Ox 
the equation of continuity [Eq. (2)] is automatically 
satisfied. 
Changing the independent variables x and y to x and 


Ou _ Ho oO (utor o) 
Ox Po oy F 


y, we have 
(14) 


By von Karman and Tsien’s transformation,” we intro- 
duce a new independent variable 


7; = (y Up) (ZL lary) ve (15) 
where L is a characteristic length. Write 
u* = u/uo (16) 


Then Eq. (14) becomes 
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Velocity distribution at the exit of a two-dimensional 
nozzle. 


Fic. 2. 

































































Fic. 3. Velocity distribution in the jet from a two-dimensional 
nozzle exhausting into a uniform stream. 
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Fic. 5. Velocity distribution in the mixing zone of two uniform 


streams. 
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Data of the Jet Flow and the Surrounding Stream 


General Data: 

Fluid in the jet and the surrounding stream 
Gas constant, R 

Specific heat at constant pressure, Cp 
Ratio of specific heats, Cp/ Cy 
Gravitational acceleration, g 

Pressure, P 


SCIENCES—AUGUST, 1949 
TABLE 1 ; ; 
Air 
53.3 ft. per ° RF. 
0.24 
1.40 


32.17 ft. per sec? 
$10 Ibs. per sq. ft. abs 


Case I: Hot Jet in Cold Stream 
Jet Surrounding Stream 
Temperature T; = 660°R. To = 330°R. 
Density pr = 0.00036 slug per cu-ft. po = 0.00072 slug per cu-ft 
Velocity uy = 1,684 ft. per sec. uw = 1,540 ft. per sec. 
Case II: Cold Jet in Hot Stream 
Jet Surrounding Stream 
Temperature T; = 330°R. T> = 660°R 
Density px = (0.00072 slug per cu.ft. po = 0.00036 slug per cu.ft 
Velocity u; = 1,684 ft. per sec. uw = 1,540 ft. per sec. 
Case III: Large Velocity Difference 
Jet Surrounding Stream 
Temperature 7, = 290°R Ty, = 330°R. 
Density pr = 0.00082 slug per cu.ft. po = 0.00072 slug per cu.ft, 
Velocity u;, = 3,000 ft. per sec. uy = 2,000 ft. per sec. 
Case IV: Case for Which the Assumption of Small Perturbations Is Valid 
Jet Surrounding Stream 
Temperature TT; = 290°R. To = 330°R. 
Density pr = 0.00082 slug per cu.ft. po = 0.00072 slug per cu.ft 
Velocity ur = 1,684 ft. per sec. uy = 1,540 ft. per sec 
n du* d Meee du* (17) For the case of mixing of two uniform flows with the 
— = u ( EY se - 
2 dy dn wp dn boundary condition of Eqs. (11) (Fig. 4), we have 


This equation can be solved by the method of succes- 
sive approximations, since p* and yu* are functions of 
u* only, as shown in Eqs. (4), (5), and (6). The 
method of successive approximations can start with 
the solution of small perturbations. 

Write 


P oe 


a” = (18) 


l + u,* 


From the solution of small perturbations, we have 


f(n); 


Substituting Eqs. (18) and (19) 
following results are obtained: 

For the case of the flow from 
nozzle under full expansion with the boundary condi- 
tion of Eqs. (9), (Fig. 2), we have 


= fi(n) (19) 


into Eq. (17), the 


oe 
u,*u*p* 


a two-dimensional 


uy (fe F ; 5 2 F in) ” 
> r ; ( ZU) 
-~ we fen” oe 74RD 
where 
ee Ee ie “ re | i 
F = exp \ [ndn/2( f + fi)] f (21) 
0 
m = (~i/uo) (L/a*x)'” (22) 
where ¥, = value of y at x = 0 and y = 1| and the 


value of ¥ is taken to be zero at x = O and y = 0. 


c= f w+ 


fi)|dn (23) 


uy l 


n F 
= d 
uy G I St+hi : 


Having computed the final u*, the y corresponding 
to u* can be calculated from the following equation: 


(24 


~ 


” 
y — Yo = (a2’x/L)' | (dn/p*u*) (25) 
0 


where y = yo when n = 0. 

To illustrate the effects of temperature and of large 
velocity difference on the properties of the jet mixing, 
some numerical calculations will be carried out here. 
The data for the various cases calculated are given in 
Table 1. 

The velocity distributions for the mixing of the flow 
from a two-dimensional nozzle under full expansion 
with the surrounding stream and for the mixing of two 
uniform streams of the various cases given in Table | 
are plotted in Figs. 6 to 10 at a typical distance down- 
stream. 

Fig. 6 shows the effects of temperature on the jet from 
a two-dimensional nozzle. For a given u0, the (1/10) 
0 decreases as the temperature of the flow in 
the jet increases. The divergence of the jet is af- 
fected little by the change of the temperature. It 1s 
probably due to the fact that, near the outside of the 
boundary of the jet, the difference of temperature and 
velocity of the jet flow from those of the surrounding 
stream are so sinall that the conditions for all the three 
cases considered are nearly the same. 
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Fic. 6. Temperature effect on the jet mixing of flow from a two- 
dimensional nozzle exhausting into a uniform stream. 


If we plot ;/(#),-0 versus y, we have Fig. 7, which 
shows that for the hot jet the rate of decrease of velocity 
in y-direction is less than that of the cold jet. 

Fig. 8 shows the temperature effects on the mixing of 
two uniform flows. The divergence of the mixing 
increases as the temperature of the jet increases. 

Fig. 9 shows the effect of velocity difference at the 
exit of the two-dimensional nozzle from the surrounding 
stream on the velocity distribution in the jet. 

Fig. 10 shows the effects of velocity difference on the 
velocity distribution in the mixing zone of two uniform 
flows. 

From the above results, we see that the properties of 
the jet mixing depend mainly on the momentum of the 
jet. The effects of the change of momentum due to the 
change in velocity and the change in temperature—i.e., 
the density—are the same. Compressibility has a 
considerable effect on the properties of the jet mixing. 

The temperature distribution can be calculated from 
the velocity distribution curve by the help of Eq. 
(4). A typical curve is shown in Fig. 11. The tem- 
perature distribution curve is a little wider than the 
corresponding velocity distribution curve. 

In all of the calculation given above, it was found 
that the second approximations give the results with 
sufficient accuracy. 


THE TURBULENT JET MIXING 


Reichardt’s Theory of Free Turbulence 
The turbulent shearing stress may be expressed as 
tT = pe(Ou/Oy) (26) 


where ¢ is the coefficient of eddy kinematic viscosity. 
According to Reichardt’s theory of free turbulence,’ 
we write 


e= R(tu. — Uuin.)O (27) 


where } is the width of the mixing zone and k is a 
dimensionless proportional factor. This theory is 
based on the assumption that the exchange over each 
cross section of the zone of mixing is constant. 

In the case of mixing of two uniform streams, Umar 


Umin, iS a constant. By dimensional analysis, 5 is pro- 


A COMPRESSIBLE 


FLUID 467 





SOLO JET 





cic 
i 
= 
_— 












10 
i\ — —NOrHEAT TRANSFER | 
—_+_ +ffs—+ ; | 
i |} ——wWOT JET 
OB 4 } 
i Br-o010 
os 
ae ee: + —{——Ag 
oa}-——+- + = 


Oo 2-—+—+>—+—_+>—_ +++ 4 


= 2 oo on ee pAb 





| 

2 Y 

Fic. 7. Temperature effect on the jet mixing of flow from a two- 
dimensional nozzle exhausting into a uniform stream. 
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Fic.9. Effect of large velocity difference on the jet mixing of 
the flow from a two-dimensional nozzle exhausting into a uni- 
form stream. 
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Fic. 10. Effect of large velocity difference on the mixing of two 


uniform streams. 
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Fic. 11. Temperature distribution in the jet of flow from a two- 
dimensional nozzle exhausting into a uniform stream. 


portional tox; hence, 


€ = €o(X L) (28) 


where € is a proportional factor whose dimension is the 
same as the coefficient of kinematic viscosity. 

In the case of mixing of a flow from a two-dimensional 
nozzle with the surrounding stream of slightly different 
mean velocity and temperature, the conservation of 
momentum in the x-direction, neglecting the higher 
order terms, gives 


b/2 
u,dy = 
—b/2 


If we assume the similarity of velocity profile at various 
sections of x, we may conclude that 


constant (29) 


€ = © (30) 
In general, the eddy kinematic viscosity coefficient 

may be written in the form 

€o(x/L) sg (3 1 ) 


eEe=> 


where 1 lies between 0 and 1 and the exact value 
of m depends on the condition of mixing. € is the 
empirical constant that can be determined from experi- 
ments. 


Equations of Motion and Solutions 


The equation of motion for two-dimensional turbu- 
Ou Ou 
ai— + w- 
x 


lent jet mixing will be 
(=) o ( o*) 
= € i — — 
0: oy rr oy ‘ oy 


Under the assumption of small perturbation, Eq. 


(32) 


(32) becomes 
Ou, x y On ‘ 
eee abel s oe 33) 
” Ox . ( oy? - 
In terms of variables x and y, Eq. (32) becomes 
ou (=) ra) ( 7 a) 
— = € _ —— 2 onan ‘ 
9 AG oy p°°u oy (34) 
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Now we introduce the new independent variable 


A BT? Fe) 35) 
for x, Eqs. (33) and (34) become, respectively, 
Ou; OX = (€)/Uo)(0?m,/dy?) 36) 
and 
Ou > ( 7 ~t) i 
—~ = & p -u oO 
OX oy oy 


Eqs. (36) and (37) are identical to Eqs. (8) and (14), 
respectively, provided that € is used for uo/po, X is used 
for x, and p* is used for y*. 

It should be noted that Eq. (4) still holds true for 
turbulent flow as shown by Frankl.’ Hence, the solu- 
tions of the laminar jet mixing can be applied to the 
problem of turbulent jet mixing of the same boundary 
conditions, provided that proper characteristic con- 
stants e) and ” are chosen. 


CONCLUSIONS 


The mean velocity distribution of a two-dimensional 
jet of a compressible fluid exhausting into a uniform 
stream are calculated for both the laminar and the 
turbulent flow. These velocity distributions are of 
similar form but with different characteristic constants, 
such as coefficient of kinematic viscosity. 

Under the assumption of small perturbation, the 
velocity distribution can be expressed in terms of 
Gauss’s error integrals. 

Exact solutions can be obtained by successive ap- 
proximations starting with the solution of small per- 
turbations. 

The divergence of the jet mixing is increased with the 
coefficient of kinematic viscosity and decreased with 
increase of free-stream velocity. The divergence in- 
creases with x-distance parabolically. 

Compressibility has a considerable effect on the 
properties of the jet. The properties of the jet depend 
mainly on the momentum of the jet. The effects of 
the change of momentum due to the change in velocity 
and the change in temperature—i.e., density—are the 
same. 

In turbulent flow, Reichardt’s theory of free turbu- 
lence has been extended to the case of compressible 
fluid. The coefficient of eddy kinematic viscosity 
should be determined experimentally. 
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SUMMARY r = wing chord, in. 
i” , ses P dei Vi = nondimensional circulation at station 7 (y, = 
The use of matrices facilitates the solution of the spanwise lift or, /bv) 
distribution for the case of a flexible wing for all conditions of LO ; - een 
: ‘ ; é ; q = dynamic pressure, lbs. per sq. in. 
load factor, incompressible air speed, aileron and flap deflection, ¥ : : 
‘ ‘ aes s aie ay = angle of attack at station 7, rad. 
and roll in a convenient form within the range of linear variation F , ah 
: : : ns , a ( ) = a square matrix 
of section lift with angle of attack. The method yields the critical a a ( 
: 3 : : , a = a matrix column 
values of wing divergence speed, aileron reversal speed, and wing- i 
aileron divergence speed, without additional computational effort, Bd =,sum of the elements of the matrix column 
for the conventional wing and may be adapted to the sweptback 0 
wing. Li = lift at station 7 (positive upwards) acting over 
interval Ay,;, lbs. 4 
, a ao = root angle of attack measured from angle for 
INTRODUCTION aes 
zero wing lift, rad. 
. roar . , “ ae , p 
HE USUAL PROCEDURE for obtaining basic flight > built in twist, rad. 
eae : f A Aa = torsional deformation, rad. a 
criteria for a flexible structure is to solve the lifting a ar . 
; : ze ‘ Cores ; P T; = chordwise torque at station 7, in.lbs. R 
line equations for an assumed spanwise distribution of é; = effective angle of attack at station i due to 
angle of attack, obtain the lift, permit the lift to twist aileron deflection of 5,[6 = (Oa/d5,)5,], rad 
the airfoil, obtain a new lift distribution for the altered pb/2v = helical angle due to roll, rad. tl 
angle of attack, and repeat the procedure until con- . = load factor in g's as 
: aay : 04; = torsional influence coefficient for twist at station al 
vergence has been obtained. This process is repeated ; ee Pop 
i 4p —_ 2 i due to a unit torque at station j, rad./Ib. in 
for each desired condition of air speed, control surface W, = weight of wing interval Ay, acting at station y, th 
deflection, load factor, etc. lbs. 
The use of matrices reduces the iteration procedure Cme.e. = wing moment coefficient about the aerodynamic 
to a single operation and permits the solution of equa- ooaner ; 
: y ; et : . My = torque due to dead weight, in.Ibs. 
tions in a convenient unit form applicable for all values : <egae 
ap : : i M; = torque due to aileron deflection, in.lbs tic 
of incompressible air speed, load factor, ete. M.n = torque due to wing camber, in. Ibs. 
The method requires certain given information, 
which may be obtained experimentally or analytically, 
as follows: 
(1) Spanwise variation of section lift curve slope. rel 
(2) Chordwise variation of aerodynamic center 
location along the span. 
(3) Spanwise variation of section moment coefficient Fo 
about the aerodynamic center. 
(4) Spanwise variation of section center of pressure | 
for a unit aileron deflection. : 
(5) Spanwise location of elastic axis. Mr 
(6) Estimation of 0a/06, for a unit aileron de- ' ! tha 
flection. ana’ “1 ' a 
These parameters are required for any of the con- x oe 
ventional methods of basic flight calculations and repre- 
sent no additional effort peculiar to this method. 
NOTATION 200 “— { 
I 
y = coordinate in the spanwise direction measured 
from the centerline, in. a, 
x = coordinate in the chordwise direction, in. 60° 
b = total wing span, in. - 
paar nae r ey . 
Received February 9, 1949. ‘ — ) a 
‘ . é F 0 
* Design Engineer. Fic. 1. Plan form and Multhopp stations. 
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(4! = column of elements equal to unity 
{1} = column of elements equal to unity or zero to | | | 
designate aileron deflection P 
(e), (m), (r) = diagonal matrices defined in the text | 
S = wing area semispan, sq.in. | 
; 
Subscripts 4 Loan Ar St. (752 
5 = symmetric S 
a = antisymmetric 
' g 
ic. = aerodynamic center é | 
. : 3 ft Loao At Sia, 141.4 
ea. = elastic axis as 
c.g. = center of gravity 8 
c.p. = center of pressure of lift due to aileron deflection a 
2 [— 
Loao Ar Sra. 100.0 
Tue Lirt DISTRIBUTION 
/ 
Let the lifting line be divided into 2m + 1 stations Loao Ar Sta. 51.8 
defined by: 
° ° — 
b wl ? ° 50 j00 /50 200 
¥; = - cos (z = 0; 1, 3, 2m) (1) 
. 2 2m 
‘ - SPAN - IN. 
as shown in Fig. 1. ae. SPAN IN den, 
Fic. 2. Twist for unit torque at wing tip. 


These are the stations at which the lift will be com- 
puted and at which the lift may be considered to act as 
aset of concentrated loads at the aerodynamic centers. 
By means of the Multhopp equations it is possible to 
obtain two matrices (following the method of reference 
1) which relate the nondimensional circulation, y;, at 
the Multhopp stations in terms of the distribution of 
angle of attack at these stations. 

For the symmetric case, involving m equations (for 


the semispan), 


(By)iYsib = (asil (2a) 
For the antisymmetric case, involving m — 1 equa- 
; g | 
tions (for the semispan), 
(Ba)t Yai aii | ais (2b) 


The lift at any station, following reference 1, is 
related to the circulation at the station by 


L; = byqAy; (3) 
For m sufficiently large, (m > 4), 
Ay; © (br/4) sin (ai, 2m) (3a) 


Introducing Eq. (3a) into Eqs. (2a) and (2b) replaces 
Mr. Benscotter’s circulation matrix (B) by a lift matrix 
that has the advantage of being completely symmetric* 
and considerably shortens the problem of inverting the 
matrix to obtain lift in terms of angle of attack. 


(B,)(bAy,) "{Lsi/g} = Lai} (2c) 
(B,)(bAy;)—"{La:/q} = {avai} (2d) 

Inverting the matrix, 
[Lsi/q} = (bAy,) (By) {asi} = (Qs){ cys} (4a) 
(La:/g) = (bAy;)(Ba)“"Laai} = (Qa)faai} (4b) 


* This was first brought to my attention by Raymond Vachss, 
of Republic Aviation Corporation. 


THE DISTRIBUTION OF ANGLE OF ATTACK 


The section angle of attack at any station along the 
semispan may be written as a linear combination of its 


components. Thus, for the symmetric case 
Qsi = Ao +t an, + Aas; (5a) 
For the antisymmetric case 
Qa = 6; + Aag; (py;/v) (5b) 


In matrix form, Eqs. (5a) and (5b) may be repre- 


sented as columns: 


{| as; = Qo} 1} + | Apt + ; Aa;;} (6a) 
jl \ pb | mt | 
{Qaii = 46 Vof + | Aaa} — Oy 7cos a f (6b) 


The elastic deformation twist at the Multhopp sta- 
tions may be expressed in terms of the forces acting at 
these stations by means of a matrix of torsional influence 
coefficients. 

{Aad = (6,73 (7) 

The influence coefficients may be obtained analytic- 
ally by any of the typical interaction methods of stress 
analysis or through test by measuring the distribution 
of wing twist for a unit torque applied to the wing tip 
as shown in Fig. 2. 

The section torque is a linear combination of its com- 


ponents; for the symmetric case 


T 31 = Mem + Lai(Xe0. — Nac.) — MWilXen. — Xeg.) (8a) 
where 
M om = Cmac.9li? AY: 
and for the antisymmetric case 
Tor @ Led%en, — Sac.) + M,, (Sb) 
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In order to adapt the deflected aileron condition to 
this analysis, the following artifice is resorted to: 


(a) 


The sum of the rigid wing two-dimensional anti- 


symmetric lift plus down-wash is assumed to act at the 
section center of pressure along the lifting line span 


— Xeg, (7) is the diagonal matrix whose ele- 
Xep.. and m is the load factor. 


SIC Xeon 
ments afe Xe — 


SOLUTION OF THE EQUATIONS 


Combining Eqs. (4), (6), (7), and (10), Eqs. (4a) and 


188 coincident with the aileron. (4b) may be rewritten as 
‘od (b) The rigid wing down-wash for the remainder of - 
dd the lifting line is taken at the aerodynamic center. [ (@) [aof 1} + {aoi}] 
a (c) The aileron moment is then defined as the mo- (Ls: { Mey \ ; 
pies ment necessary to transfer the rigid wing antisymmetric {es + (2)(04)| a rim n(m)} Vi 4 (Ila) 
lift in the region of the aileron from the section center 

353 | of pressure to the section aerodynamic center. Lt (25) (013) (e){Lsi} al 
cmeill (d) The incremental lift due to wing twist resulting r 1 Pb ri\ 

from the application of load and moment is carried (a) {3} - }cos om t) 

entirely at the aerodynamic center. +t = 2 - (11b) 

qf + (Qq)(01;)[(r){ Lat crigiay ] 
My, = Lai (rigiay(Xac. — Xe.p.) (9) | + (@,)(0,) (Cf Lar ] 








Eqs. (Sa) and (Sb) may be written in matrix form 
. Since the only unknowns are {|L,,j and {La}, Eqs. 


: - {Mem \ ; (lla) and (1lb) may be solved by iteration. Let 
’ (Ts = q ) q ‘tes (e)iLs3 — nmi Wit (10a) {7./9}1 be the first approximation to the solution of 
: Eq. (lla), obtained by neglecting the last term on the 

{Toad = (e){Laid + (r){ Las crigiay} (10b) right-hand side of the equation. Then the incremental 


lift due to the twist resulting from the application of 
qiL./q}' at the aerodynamic centers along the wing 
span will be 


Where (e) is the diagonal matrix whose elements are 
“ea — Xa, (m) is the diagonal matrix whose elements 
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{ AL;/g}! = (2)(6;;)(e)gt Li/q}! (12) 
The second incremental lift will then be 

{ AL: /g}? = (2)(6:;)(e)gt AL: /q}! (13) 


etc. A solution exists if the series formed from the 
incremental lifts converges. Normally this would 
prove a cumbersome computational task; however, in 
matrix form a simplification results from the character- 
istic property of matrices. 

Let the matrix (2)(6;;)(e) be defined as (F). Then, 
as shown in reference 2, the product of (F){t} yields 
another column {/}!. Repeating the operation a suf- 
ficient number of times transforms the column into 
the dominant “eigen function” of the matrix (/) which 
has the property that 
(14) 
where \ is the dominant “eigen value’’ of the matrix, 
(F). 

As a consequence of Eq. (14), 

[> fF) NA = [SS ‘A = 1/0 — a ]f4 (14a) 
k=0 k=0 

Once the arbitrary column | AL; q\" has been trans- 
formed into a “eigen function” of the matrix (Q)(6;;)(e), 
the sum of the series of the incremental lifts may be 


written in closed form 


= TW ALi? Aetee . 
¥T(Q)(0,;)(e) | > = —. (15) 
| 59 (2) (64)(0)] lq q l—rdA\qlqJ 


where /: is the number of iterations required for the 
ratio of the corresponding elements in succeeding col- 
umns to yield the proper constant, X. 

The physical significance of \ is immediately appar- 
ent, since the matrix g(2)(6;;)(e){1} yields the incre- 
mental lift distribution for unit loads applied along the 
wing aerodynamic centers. 

Wing divergence will occur when an applied unit load 
distribution along the aerodynamic centers produces 
an incremental load distribution equal to the applied 
load in both sign and magnitude. Wing divergence is 
given by 

1 — gars = O 


12V'391 A, m.p.h. (1.A.S.) (16) 


Similarly, the value of \, corresponds to the anti- 
symmetric wing-aileron divergence dynamic _pres- 
sure. 

For the case of the Multhopp equations, the itera- 
tion process required to bring the arbitrary matrix 
column to the canonical form is of the order of one or 
two iterations. This is amply illustrated in the example 
included in the paper. 

The solution of Eqs. (lla) and (11b) for the case of 


the most rapid convergence (i = 1) is given by 
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the symmetric distribution of @ and one for the 
antisymmetric). However, care must be exercised to 
obtain the proper flexibility matrix of torsional in- 
fluence coefficients. For the sweptback wing the tor- 
sional flexibility is a function of the torsional rigidity, 
the bending rigidity, and the angle of sweepback. 


For the case of a sweptback wing, the Multhopp 
equations of reference 1 are replaced by a similar lift 
distribution” matrix obtained through the Weissinger 
The Weiss- 


inger equations yield, as before, two matrices (one for 


lifting line theory as given in reference 3. 





























> 
7) Basic FLIGHT CRITERIA 
Eqs. (17a) and (17b) are the unit solutions of the lift distribution referred to above and form the basis for several 
of the important basic flight criteria. 
NNEL) From Eq. (17a): The wing C, for the flexible wing is given by 
T . fLsi\ , q = eal | 7 
ra r + ot > — 
ib aq ln, 1-mTL @ Ja 
ue. a m ‘7 1 q m TAL. 1 
VERGENCE m + > Ped J AL, \ 
pp» 0 lq ) a 1 — Asq “O q }% (18) 
. 0 ;, 8) 
¢ = = m pl 
. gS S| + q {4%} 
1 — Aq 9 qd f Mem 
n m AL; 1 
? ta p ‘ t 
a I 7 0 q J Mw : 
sible). 
" by ! implies the summation over the m values of | 
2 ( 
mm) @ The wing lift slope is given by 
95 ) 
6 966 L > {Aza 
" 1. Ss 4a “~ JAL,;t! 1 —  X, 
a ”_ Sp J F l + q x ‘ il | —_ ~ — q 07 # Mo (19) 
SLO \ q Pin 1 — Aq q fa, G.W. l f ALi! 
- 4 < { 
2 l~-aAg yd ] q Mw, 
The second factor arises from the fact that in flight the application of lift to the wing accelerates the airplane, 
causing an incremental increase in wing twist due to the moment of the wing dead weight about the elastic axis. 
This additional effect may result in a sizeable discrepancy in the estimation in C,, for an elastic wing measured in 
the wind tunnel as compared with the C,, experienced in flight. 
The induced section drag may be obtained directly from the section lift distribution as shown in reference 4. 
From the Prandtl equation 
dL(y)/dy = C,c(y)qta(y) — [w(y)/V]} (20) 
5 and 
dD(y)/dL(y) = w(y)/V (21) 
For the Multhopp intervals [ Ay; > (2b/4m) sin (wt/2m)], the section induced drag is given by 
Dg = (ao + a; + Aa — [brLsi/C,,Cigdm sin (4t/2m)]} Ls, (22) 
1 > " - one ° ° ° ° ° 
wh From Eq. (17b): The rolling moment for the semispan is given by 
Mu \ = 
7 gé veg + ; Vi 
c OLE Go tron Oa 
RM. <= 6 SAL, |! 
= = > Ladi = he q J ni vip (23) 
2 0 1 — aq “OD l qd M5 
. b m { 7 1 m A A 1 
a 2v 0 \ q ds 1 — Aad 0 l q f 2n 
1 
| The roll effectiveness (for the pb/2v per unit aileron deflection corresponding to R.M./2 = 0) is given by 
: (170) } n = [(pb/2v)/6 (elastic)]/[(pb/2v)/6 (rigid) ] (24) 
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where 
- if, . q ™ JAL,, : =. (AL i 
Ba ot + LE Sod te | 
pb : og Js — Lo q ds 0 q M5 seta 
i 5 (elastic) = - cae (25) 
2uU ie 
\— vibe Sai peer ope fom vaton 
0 (gq ) 20 —dg ol g 2r 
4 | a fle V'/ fle . \' ‘ 
—/§ (rigid) = — VW —= Var (25b) 
2y x  -_ f, x * fos 
The steady state pb/2v for a given aileron deflection 6, is given by 
pb/2v = n[(pb/2v)é (rigid) ]6,(0a/06,) (26) 
The aileron reversal speed is given by 7 = 0 
me fALa Ui fala . | 
pe nt a > i the. 
/ 0 qg 4 0 q Mé i 
Greversal = 1 Na ee a a ™ (27 
p Les . 
st ys 
0 lq F 
The induced section drag for the antisymmetric case is given by 
y bxrL, 
Du = E + Aaa a — at 28) 
v O,C:g4r sin (2;/2m) 
EXAMPLE aileron c.p. = 0.50C, 
ee , = — (2/57.3) cos (ri/12 
The method is illustrated for m = 6 for a typical — iat 8 (w4/12) 
wing having the following parameters: Siinittiiiaiatnts 
b = 400 in. ' Benscotter, S. U., Matrix Development of Multhopp’s Equa- | 
Cs = (100 — 40) cos (4/12) in. tions for Spanwise Air-Load Distribution, Journal of the Aero- 
C, = 5.9 nautical Sciences, Vol. 15, No. 2, pp. 113-120, February, 1948 
; - fm ie 2? Frazer, Duncan, and Collar, Elementary Matrices; The 
am, = —0.02 MacMillan Company, New York, 1947. 
Nec, = 0.25C, 3’ Weissinger, J., The Lift Distribution of Swept-Back Wings, 
cars 0.40C N.A.C.A. T.M. No. 1120, Translation, March, 1947. 
Xea, = 0.400, * Von Mises, R., Theory of Flight, pp. 231-244; McGraw-Hill 
Xex, = 0.50C, Book Company, Inc., New York 
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A Formulation of the Aeroelastic Problem 
for a Swept Wing’ 


JOHN W. MILESt 


University of California at Los Angeles 


SUMMARY 


On the basis of certain simplifying structural and aerodynamic 
assumptions, notably that wing bending and torsion occur along 
a straight line, and strip theory, the differential equations of 
equilibrium of a swept elastic wing subjected to aerodynamic and 
inertial loads are formulated. Integral forms of these equations, 
representing energy balances, are also developed. Formulas for 
the lift and moment coefficients in terms of the initial configura- 
tion and elastic deflections of the wing are given. Various 
methods of solving the equations are discussed, and an approxi- 
mate solution to the energy equations is included. An appendix 
is given on the calculation of effective angle of attack distributions 
over a rigid wing 


NorTaTION 


A = amplitude of torsional deflection, cf. Eq. (38a) 

B = amplitude of bending deflection, cf. Eq. (38b) 

Gas = cf. Eqs. (47), (48) 

Cii = cf. Eqs. (49)-(52) 

CL = lift coefficient, cf. Eq. (15) 

Cr° = lift coefficient for rigid wing, cf. Eq. (20) 

Cr* = lift coefficient corresponding to a section lift coeffi- 


cient distribution 7, cf. Eq. (18), where i = f’, g, 
y, or ¢, and unit lift curve slope 
Cu; = ef. Eq. (74) 


Cu = pitching moment coefficient, cf. Eq. (16), positive 
stalling 

Cu,,, = cf. Eq. (75) 

Cy® = pitching moment coefficient for rigid wing, cf. Eq. 
(21) 

Cy' = moment coefficient corresponding to i, cf. Eq. (19), 
where 7 = f’, g, ¥, or ¢, and unit lift curve slope 

Ci = rolling moment coefficient, cf. Eq. (17), positive 
when tending to depress starboard wing 

Ci = rolling moment coefficient for rigid wing, cf. Eq. (22) 

E = flexural modulus of wing 

Ex = cf. Eq. (46) 

F = vertical load on wing, cf. Eq. (5a) 

G = torsional modulus of wing 

Gu = cf. Eq. (45) 

I = structurally effective moment of inertia of wing 


section (perpendicular to elastic axis) with respect 
to neutral axis (chordwise) 


bY = structurally effective polar moment of inertia of 
wing section with respect to elastic axis 

ZL = total lift on wing 

= surface area of wing 

T = twisting moment about elastic axis, cf. Eq. (5b) 
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ct 
c)° 


Cy 
Cm 


c.g. 
c.p. 


ki, My 


A.R. 
l.a.c. 
m.a.c. 
m 


free-stream velocity—i.e, airplane’s steady flight 
velocity along x axis 

distance of airplane center of gravity forward of the 
aerodynamic center of the mean aerodynamic 
chord, the latter being defined for the rigid wing 


= wing span, measured tip to tip 


chord, measured parallel to plane of symmetry 

chord, measured perpendicular to e.a., cf. Eq. (3) 

section lift coefficient, cf. Eq. (6a) 

section lift coefficient for wing in absence of aero- 
elastic effects 

section lift coefficient due to source 7, cf. Eqs. (26) 
and (27) and the Appendix 

moment coefficient for section perpendicular to 
l.a.c., cf. Eq. (6b); also m.a.c. 

center of gravity 

center of pressure 

spanwise (y) coordinate of inboard end of control 
surface 

eccentricity—i.e., fraction of chord of e.a. aft of 
l.a.c. 

elastic axis 

normalized bending deflection, cf. Eq. (37a) 

normalized torsional deflection, cf. Eq. (37b); also 
acceleration due to gravity 

bending deflection of e.a., positive up 

subscripts 

radius of gyration of wing section (perpendicular to 
e.a.) with respect to section c.g. 

ef. Eqs. (71)-(73) 

length of e.a. over semispan [{/ = ('/2)b sec y] 

aspect ratio 

line of aerodynamic centers 

mean aerodynamic chord 

mean lift curve slope, corrected for compressibility, 
sweepback, and induction effects, cf. Appendix 

angular velocity in roll about x axis, positive for star- 
board wing going down 

angular velocity in pitching about y axis, positive 
for nose going up (stalling); also dynamic pressure 


= g/qmaz. (telative dynamic pressure) 


angular velocity in yawing about z axis, positive 
for starboard wing going aft 

coordinate measured along e.a. from 0 at root to 1 
at tip 

time 

perturbation in forward (x axis) velocity of c.g. 

sideslip (y axis) velocity of c.g. 

plunging (z axis) velocity of c.g. 

axis fixed in direction of flight, origin at c.g. 

distance of a.c. of m.a.c. forward of intersection of 
the l.a.c. with the root chord, generally negative, 
cf. Fig. 1 

a+xXxm+ ytan y 

axis perpendicular to plane of symmetry, positive to 
starboard, origin at c.g. 

axis perpendicular to x and y axes, positive down, 
origin at c.g. 
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A = cf. Eq. (55) 
I = cf. Eq. (61) 
2 = cf. Eq. (64) 


= angle of attack (incidence) of section parallel to 


Qa 
plane of symmetry 

a = angle of attack distribution due to source 7, cf. 
Eqs. (28) to (35) and the Appendix 

B = dihedral angle, measured as a rotation about x axis 
(i.e., viewed along y axis); the complement of the 
dihedral angle between the plane of symmetry 
and the tangent plane to the wing at a given sec- 
tion 

¥ = the sweepback angle of the l.a.c. with respect to the 


y axis; also the sweepback angle of the e.a. on 
the basis of the preliminary assumptions 
5 = angle of rotation of control surface about hinge line, 
positive for starboard surface deflected down 
€ = cf. Eq. (65) 
= angular rotation of wing section perpendicular to 
e.a. due to motion of air frame as a rigid body 


= distance of section c.g. ahead of e.a., measured per- 


"eo 


pendicular to e.a. 
0 = angle of pitch of airplane about y axis 


cf. Eq. (42) 


K = 

» = cf. Eq. (41) 

rm = cf. Eq. (43) 

v = cf. Eq. (44) 

p = mass density of air 

a = weight per unit (swept) span 

’ = angle of roll of airplane about x axis 

¢ = angle of twist of section perpendicular to e.a. 

¥ = angle of yaw about z axis; also bending slope 
(Oh/Ods) 

g = dimensionless span coordinate, cf. Eq. (4) 

( )? = ( ) evaluated for rigid wing 

( )o = ( ) evaluated at root section 

®. = d( )/dt 

t.}’ = derivative with respect to argument of ( ) 

( )e = O( )/ds 


INTRODUCTION 


5 inn ADVENT OF THE SWEPTBACK WING has changed 
considerably the aspect of the aeroelastic problem 

i.e., the calculation of the aerodynamic load distribu- 
tion over an elastic wing. The primary factor in the 
analysis of a swept wing will usually be the wind in- 
cidence introduced by bending, which affects only the 
dihedral of an unswept wing. In the case of an unswept 
wing, torsional deflection generally shifts the loading 
outboard, with a resultant tendency toward wing diver- 
gence, although control surface loading (being aft of the 
elastic axis for normal structures) produces deflections 
that decrease the outboard wing loading, with a result- 
ant tendency toward control surface reversal. On the 
other hand, bending due to either normal or control 
surface loading of a sweptback wing tends to shift the 
aerodynamic loading inboard; conversely for a swept- 
forward wing. Accordingly, the possibility of wing 
divergence for a sweptback wing is decreased, while the 
possibility of control surface reversal is increased; 
conversely for sweepforward. Moreover, elastic de- 
flection of swept wings will produce severe c.p. shifts, 
tending to produce low (or negative) static stability for 
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a sweptback wing, excessive static stability for a swept- 
forward wing, and undesirable trim requirements in 
both cases. 

The dynamic stability is also affected by aeroelastic 
effects because of the changes in the aerodynamic 
derivatives. The following analysis will be applicable 
to the calculation of these stability derivatives if it is 
assumed that the natural frequencies of vibration of the 
elastic airplane are large compared with the natural 
frequencies of the motions considered in dynamic sta- 
bility and that aerodynamic lag may be neglected. 

The most difficult part of the aeroelastic problem is 
the prediction of the elastic deflection of a swept wing 
under a specified loading. While a considerable amount 
of experimental work should be carried out in order to 
facilitate rational assumptions, it appears that the 
rather drastically simplifying assumptions made below 
may suffice for wings of large panel aspect ratio (A.R. 
sec* y), at least for preliminary design purposes. (It 
may be remarked that the full aerodynamic benefits of 
sweepback cannot generally be realized without going 
to reasonably large aspect ratios.) However, the 
analysis to be given would have to be extensively modi- 
fied before it could be applied to such configurations as 
the delta wing. Fortunately, it appears that the 
aeroelastic problem will not be too serious for the delta 
wing. 

The aerodynamic portion of the aeroelastic problem 
also presents difficulties, but they are similar to those 
encountered in the aerodynamic analysis of a rigid 
swept wing and may be handled in the same manner. 

The following analysis formulates the basic equations 
for the calculation of the lift distribution on an elastic 
wing and the associated lift and moment coefficients. 
Having these equations, the aeroelastic analyst may use 
whatever methods of solution he deems appropriate, 
and certain special solutions have been carried out in 
the papers listed in the references. A general survey 
of the field of aeroelasticity, together with an extensive 
bibliography, has recently been given by Collar.' 


ASSUMPTIONS 


The following assumptions will be made throughout 
the analysis: 

(1) The aerodynamic centers of the chord sections 
lie on a straight (spanwise) line (l.a.c.), and these chord 
sections will be defined parallel to the air-frame plane 
of symmetry. 

(2) There exists an “elastic axis” (e.a.), and all sec- 
tions perpendicular to the e.a. remain parallel during 
A torsional stiffness distribution 


torsional deformation. 
The angle between the 


GJ, is assumed along the e.a. 
l.a.c. and e.a. is assumed negligible. 

(3) The bending of the wing can be represented by 
the bending of the e.a. A bending stiffness distribution, 
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wept- : . ° ° . 
te 3 EI, calculated for sections perpendicular to the e.a., x 
its in 
: is assumed along the e.a.* x 
H (4) Chordwise deformation is neglected. In partic- C6 
lastic § ular, the control surfaces are assumed to deflect (elas- -y +—-Z 
lamic tically) as an integral part of the wing structure, so that 1 
cable no section moment coefficient (with respect to the l.a.c.) 
it 1s is introduced because of elastic deformation. a 
ot the (5) The geometrical angle of incidence (a), for a 
tural chord section at a given span station (£), due to elastic 
> sta- deflection is given by 
| ae 
a(é) = ¢(£) cos y — [0A(E)/Os] sin y (1) 
“Ill 1S 
wing where ¢(£) is the angle of twist of a section perpendic- t | 
ount ular to the e.a. at &, Oh(E)/Os is the bending slope of ecc 
er to the e.a. measured along the e.a., and y is the sweepback | { 
the angle of the wing [identical for e.a. and l.a.c., due to ie \ | 
elow (1) and (2) above]. No allowance is made for the fact Ss hk 
ALR. that chord sections taken parallel to the air-frame plane Secrion A-A 
(It of symmetry (in which a is measured) and perpendicular eCTI 
its of to the e.a. (in which ¢ is measured), at the same station ' 
roing fon the e.a., intersect the l.a.c. at different stations. 
the (6) The geometrical dihedral (8) due to elastic B 
10di- deformation is given by ~y 
ives B(é) = o(£) sin y + [0h(E)/Os] cos (2) | 
the Z 
It: This dihedral is measured as the slope of the wing alon : é , 
lelta : i ‘ I 8 § Fic. 1. Three views of a swept wing. 
a line perpendicular to the air-frame plane of sym- 
metry. . . 
blem cs Kids : ELASTIC EQUATIONS 
; (7) The chord c’(&) measured perpendicular to the 
10se , ; ' . , ian a 
‘oid e.a. at a given span station, &, is related to the chord Let (OF /0s) be the applied force (positive in the 
rigic A i : : see Aggie . : 
S c(t) defined in assumption (1) by lifting direction) per unit length of the e.a.; (O07'/Os), 
er. . ne - ae 
: the applied twisting moment (about the e.a., positive 
ions c’(£) = c(&) cos y (Ss) . ; : : ‘ 
Aon: in the stalling sense); /, the bending deflection of the 
astic (S) The effects of drag are neglected. e.a. (positive up, 1.e., along the negative z axis); y, the 
nts. twist about the e.a.; and EJ and GJ, the flexural and 
howl CooRDINATES torsional stiffnesses in a section perpendicular to the 
ate, e.a. Inertia loads, if present, are assumed to be in- 
it in The coordinates to be used are depicted in Fig. 1. eJyded in F and T by virtue of D'’Alembert’s principle. 
rvey x is measured in the direction of flight, y is measured The partial differential equations of elastic equilibrium 
sive perpendicular to the plane of symmetry and positive to _ petween the internal stresses and external loads are then 
starboard, z is measured positive down, and s is meas- given by 
ured along the e.a. x, y, are stability wind axes (i.e., : 
directions fixed in space) with their origin at the c.g., (0*/Os*) [ET(0*h/ds*)] — (OF/0s) = 0 (va) 
whereas the s axis is a body axis with its origin at the (0/0s)[GI(de/ds)] + (O7T/ds) = 0 (5b) 
lout intersection of the e.a. with the plane of symmetry. ; ; ; 
It will be convenient to define the dimensionless span Ina strict sense, Eqs. (5a)and (5b) are integrodifferential 
—_ coordinate equations, since the applied aerodynamic loads at a 
ail given station depend on the geometrical incidence over 
y) M9en - . > 
— e_ (27) . (27) _ fs the entire span. 
ane g= = = (4) a ‘ ae : 
b b l rhere will be two principal sources of loading, aero- 
/ : , . dynamic and inertial, where gravitational loading may 
en and / = (1/2)b sec y, where 0 is the tip-to-tip span and a ' 5 ‘ ie may 
: iieeeiete _ : ‘ , be included in the latter category. It is convenient to 
ring represents the fraction of the semispan from the root aw 6 load; : le 
a ; : consider these loadings separately. 
rion chord to the station y and is equal to unity at the star- poem ‘ 
the board tip. £ has the virtue that it refers equally well 
10 the « ; : AERODYNAMIC LOADS 
0 the swept (s) or nonswept (y) spanwise axes. 
: , a ; et c, be the section lift coefficient < a > sectio 
by * Assumptions (2) and (3) appear to be supported by the ex- Let c, be the : ection lift coefficient and ¢ 59 the sect . 
on, perimental results of Zender and Libove, NACA T.N. No. 1525, ™oment coefficient (about the L.a.c., positive stalling). 
1948, That part of c; which is due to the mean incidence of a 
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section is presumed independent of whether the section 
is parallel to the plane of symmetry or perpendicular to 
the e.a., because of the basic geometrical assumptions. 
On the other hand, c,, and that portion of c,; which de- 
pends on camber are intrinsically dependent on the 
orientation of the section. Inasmuch as the most im- 
portant source of camber contemplated in aeroelastic 
studies arises from a rotation of a control surface about 
a swept axis (assumed parallel to the l.a.c. and e.a.), 
camber will be measured in sections perpendicular to the 
e.a. On the basis of the last assumption, the applied 
aerodynamic loads are given by 

OF/Os = qc’ c, = qe Cc, cos 7 (6a) 
OT Os = ge’?(Cm + eC) = Ge?(Cm + eC;) cos? y (6b) 


where g is the dynamic pressure (pU/*/2) and e (the ec- 
centricity) is the fraction of the chord of the e.a. aft 
of the lL.a.c. 

Because of the assumption of no chordwise deforma- 
tion, c,, is unaffected by ¢ and h, while c, is modified 
by the change in geometrical incidence. If a represents 
the distribution of geometrical incidence in the absence 
of elastic deformation, then the addition of Eq. (1) and 
the introduction of Y = Oh/Os gives the total geo- 


metrical incidence as 
a = a + gcosy — ysin y (7) 


Corresponding to each of the three terms on the right 
side of Eq. (7) is a lift distribution, and the total lift 
coefficient distribution may be written 


) p ° > 
c, = c} + cf cos y — cf sin (8) 


where c may now be assumed to include the effects of 
camber, as well as incidence. 

Because of induction effects, c,; at any station de- 
pends on the incidence at all stations, and the explicit 
deduction of Eq. (8) from Eq. (7) therefore depends on 
the solution of an integral equation. (Since this integral 
equation is presumed linear, the one-to-one correspond- 
ence between c! and a’, cf and yg, and cf and y is 
nevertheless proper.) While it is entirely possible to 
carry out a solution of Eqs. (5a) and (5b) in conjunction 
with this integral equation, the required labor is rather 
excessive in view of the accuracy justified by the various 
geometrical and structural assumptions on which the 
present analysis is based. Accordingly, Eq. (8) will be 
solved by strip theory, wherein it is simply assumed that 
a mean lift curve slope, m, exists such that 


c, = + m(g cos y — ysin y) (9) 


The most expedient choice of m is probably the mean 
lift curve slope (OC,°/O0a°), as determined from calcu- 
lations or wind-tunnel measurements for the rigid wing. 
However, it should be realized that m actually refers to 
a twisted wing and may differ from the value calculated 
or measured for an untwisted wing. Since the only way 
of accurately evaluating m for the twisted wing is a 
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solution of the integral equation relating Eqs. (7) and 
(8), it must generally be taken as (OC,°/0a°), at least 
for symmetrical lift distributions. In the case of asym. 
metrical lift distributions, induction effects are stronger, 
and the wing may be approximately treated as two sepa- 
rate panels, so that the effective aspect ratio is half the 
true aspect ratio. The calculation of m is discussed 
further in the Appendix. 


INERTIAL Loaps 


The acceleration of the air frame will generally pro- 
duce inertial loads over the structure, which in turn will 
produce deflections. It will be assumed that these de- 
flections are in phase-with the accelerations that produce 
them, which is tantamount to the assumption that the 
time constants or periods of the air-frame motions 
considered are large compared to the natural periods of 
the wing vibrating as an elastic system. On the basis 
of this assumption, it is not necessary to account for 
the inertial loads corresponding to or ¢. (Similarly, 
the direct aerodynamic loads due to h and ¢ are neg- 
lected.) In those cases where such assumptions are 
not justified, there is essentially coupling between the 
motions of the airplane generally included in the dy 
namic stability and flutter situations, respectively, and 
the aeroelastic problem becomes rather considerably 
more complex. The study of such problems will there 
fore be delayed until a later date. 

Let o be the weight per unit span (along the swept 
axis) of the wing, let ¢ be the distance of the center of 
gravity (c.g.) of a given chord section (perpendicular 
to the e.a.) forward of the e.a., and let & be the radius 
of gyration of the same section with respect to its c.g. 
The net vertical acceleration of the section c.g. will be 
taken as 2 (positive along the z axis), while the net an- 
gular acceleration of the section will be taken as jj, 
positive in a stalling sense. 

It should be specifically pointed out that Z is not 
identical with w, the z component of the airplane c.g. 
acceleration. Moreover, 2 includes a term —g, corre- 
sponding to normal static loading of the wing due to its 
own weight. 

In terms of the parameters defined in the last para- 
graph, the inertia loads are given by D’Alembert’s 
principle as 

OF Os = (10a) 


OT /ds = 


(o/g)2 


(ao g)(t% — Rk?) (10b) 


EXPLICIT DIFFERENTIAL EQUATIONS 


The differential equations can now be explicitly stated 
in terms of the dependent variables by substituting the 
aerodynamic and inertial loads given by Eqs. (6) and 
(10), respectively, together with the strip theory lift 
coefficient of Eq. (9), in the basic equations [Eqs. (5a) 


and (5b)]. Thus, 
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(0°/ds*) [EI (Oy /Os)] + gm sin y cos y cy — 
2 p _.0 on 
gm cos’ ycg = q cos yer; + o(%/g) (Ila) 


—gme sin y cos® yc*°y + (0/0s)[GJ(0¢g/0s)] + 
qme cos® yc*¢ = —g cos* yc?(ec,® + Cm) + 
(a/g)(k? — £2) (11b) 


Eqs. (lla) and (11b) are simply a statement of equilib- 
rium for a section of the wing c’ ds(= c cos y ds). 
The flexural equation [Eq. (lla)] gives the required 
balance between the flexural resistance force —0[(EJ 
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incidence produced by bending, the aerodynamic force 
gmc’ ds ¢ due to incidence produced by twisting, the 
applied aerodynamic load gc’ ds c,°, and the inertia load 
(o/g) ds on the element c’ ds. Similarly, the torsional 
equation [Eq. (1lb)] gives the required balance be- 
tween the applied moments about the elastic axis, and 
the terms, in order, represent the aerodynamic twisting 
moment due to torsional strain, the twisting moment 
due to incidence produced by twisting, the twisting mo- 














ment due to the applied aerodynamic loads, and the 


h,,),], the aerodynamic force —gmc’ ds ¥ due to _ inertial twisting moment. 


ENERGY EQUATIONS 


It is frequently convenient to deal with an integral form of Eqs. (lla) and (11b), representing energy equations 
for the wing. Thus, if both sides of Eq. (1la) are multiplied by # ds and integrated over one wing (s = 0 to /) 
and, similarly, Eq. (11b) multiplied by ¢ ds and integrated, the stiffness integrals are integrated by parts, viz: 


2? oh ra) oh Oh __ dh! J *1( on) ; 
EI — )h ds = h—( EI — — EIx- EI ds (12a) 
f <( Os ) 7 = 45 ( on) | ds 32, 


l \l 1 9 
. 7 OY | f , oe 
I = 7J | _ J nate Ss { 12b } 
f eo(as$ yas wat 72 .* (= ” 


and the cantilever boundary conditions 





oh 0/ oh oh 

| ms = ={ er “aia «at (13a) 
: a ee = =) = | 

¢ls <0 = GJQ¢g/%s)|, -; = 0 (13b) 


assumed, the resulting equations may be written (reintroducing 0h/0s = y) 


*| i i | 
{ EI(0*h/Os*)*ds + gm sin y cos y [ ch(Oh/Os)ds — qm cos? y [ cgh ds = q cos y / ccjh ds + 
0 0 0 0 


‘| 
| a(z/g)h ds (14a) 
0 


a Oh ae | 
qme sin y cos? fe ec? — gds+ ‘as (2 4) ds — qme cos® vf c’g’ds = q cos? y / c(ec;® + cme ds + 
0 Os 0 0 0 


*| “1 
4 ok? ¢ ods -{f of (: ‘Ve ds (14b) 
0 g 0 


Eqs. (14a) and (14b) represent an energy balance for the wing. The terms on the right represent the work done 
by the applied forces in deflecting the wing. The elastically stored energies are given by 


1 1 
(1 af EI(07h/Os*)*ds and (1 » GJ(O¢/Os)*ds 
0 0 


Hence, the energy represented by the remaining integrals on the left, plus an amount of energy exactly equal to. 
the elastically stored energy (since the integrals do not contain the '/: factor), is delivered to the airstream during 
the process of deflection.* It is evident that positive bending deflection (for positive y or sweepback) tends to 
return energy to the airstream, whereas torsional deflection tends to absorb it. 

That value of g for which the energy balance on the left vanishes identically—.e., 
elastic energy is just furnished by the airstream in the absence of applied forces—corresponds to the wing diver- 
gence speed. Evidently the possibility of divergence depends on the relative amounts of bending and twisting 
and on the swegpback. Thus, a wing with high sweepback and a reasonable ratio of torsional to flexural stiff- 


_ 


the point at which the stored 


* The manner in which this transfer takes place cannot, of course, be studied through the steady-state equations, being essentially- 
a transient phenomenon. 
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ness will not generally possess a real divergence speed, whereas a wing with sweepforward and/ora low ratio of 
torsional to flexural stiffness will generally exhibit a relatively low divergence speed. 


LIFT AND MOMENT COEFFICIENTS 


The lift coefficient distribution over the wing is given by Eq. (9). 


The corresponding overall lift and moment 


coefficients (Cy, about the y axis, positive stalling, and C, about the x axis, positive when tending to depress the 


starboard wing) are given by 


ee lift i. , E . ( ny] , - 
= = ¥ os ‘OS 9— § U { }.)) 
L gS s J, Ci Ci + m\ Cos y¢ sin 135 ay ) 


Cc, am 


M, = 4% 2 | ? 
qgSb ~ Sse q «i Se m( cos ye — sin 15 


uM <-2 f"* _ oh 
a : cl c,’ + m\ cos ye — sin y fa + x, + y tan y]dy 16) 
Cn Cuts Os “ . 


oh 3 
y dy 17) 


where a is the distance of the c.g. ahead of the aerodynamic center of the mean aerodynamic chord, where the 


latter is located at y = 
l.a.c. with the root chord—cef., Fig. 1. 


It is convenient to introduce the following non- 


dimensional coefficients: 


"b/2 
c,' = (2 s) | cy dy (1Sa) 
*b/2 
Cc. = s) | cy dy (18b) 
*b/2 
Cy” = (2 os) | cy ydy (19a) 
0 
*b/2 
Cu” = (2, os) | cyy dy (19b) 
0 
*b/2 
Cc,” = 2 sf cc (dy (20) 
— b/2 
Cy" = ah cc’la +x, + ytan yldy (21)* 
Cm 0 
—_9 fe? 
C) = 2 | ccry dy (22) 
eee sooo 


In terms of these coefficients, the results of Eqs. (15) 


(17) may be written 


Cr = C,° + m(C_% cos y — C," sin y) (23) 
Cu = Ca? — (‘ t **) m (C,* cos y — Ci =>) - 
Cm 
(° — mic yw’ cos ¥ Cy’ sin y) (24) 
m 
C, = Cf — m(Cy* cos y — Cy" sin 7) (25) 


It should be remarked that the value of m used in 
Eq. (25) will generally be smaller than the value used 
in Eqs. (23) and (24) because of the asymmetrical 
induction effects. 


* For the usual definition of ¢m, Cvq, = —(a/cm)Cia,, 


—X, cot y, and x,, is the distance of the a.c. of the m.a.c. forward of the intersection of the 


INITIAL LIFT DISTRIBUTIONS 


It will generally be expedient to divide the initial 
lift coefficient into a linear superposition of its con- 
stituents. For most analyses it will be sufficient to 
include the following terms: 


u v 
c, _ e.° a Cig,” ato oe €15°6 + Eg (| ) a ro (7.) + 


w »b i rb 
Ci,” (*) + Ci” (5) > Cig (: ) +c,° (= : (26)T 


where c;, is the result of camber and initial (built in) 
twist with respect to the root section, ap is the angle of 
attack of the root section, 6 is the control surface de- 
flection (either symmetrical or asymmetrical), ™ is the 
perturbation (on Ll’) in forward (%) velocity, v is the 
sideslip (¥) velocity, w is the plunging (3) velocity, Pp 
is the rolling (¢) angular velocity, q is the pitching (6) 
angular velocity, and 7 is the yawing () angular veloc- 
ity. The positive directions of the angular velocities 
with respect to their axes are specified by the right- 
hand rule. 

If strip theory is utilized, the various section lift 
coefficients may be written 


Cc,” = ma; 


(27) 


where a; is the angle of attack distribution due to the 
source t. If induction effects are to be exactly accounted 
for, ¢,, must be obtained from a, through the solution 
of an integral equation rather than simply by multiply- 
ing by m. The various a, are derived in the Appendix 
and may be written: 


(28) 


a, = 0a/Oay = 1 


t In some cases it may be desirable to include c;,,° and ¢;,;", 
because of wing down-wash at the tail. 
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Oa ( 2) rrre 
a= == — OSs ) = 
¥5 ds a6 Jc, cos y_ (¢ 7 £ 2 
=0 (0Zy<d) 
(29) 
a, = 0a/O(u/U) = 2a (30) 
a, = 0a/d(v/U) = B+ a tan y (31) 
Ay = 0a/0(w/U) = 1 (32) 
a, = 0a/0(pb/2U) = & — sin ycos y(c/b) (33) 


Oa a a ( = =) 4 l — 
= - = an - cos? 
” O(gem/U) Cm ? Ce 2 . Cu 
(34) 
Oa _ | “y+ ‘ ( tI 
a O(rb/2U) a 


r 
[tan ya + B] — 2&a — (<) cos? y(a tany — 8) (35) 


where d is the coordinate of the inboard end of the 
control surface; (—0a/06)c; is the amount (Oa) by 
which the angle of attack of a section must be decreased 
for a positive (down) control surface deflection (06) 
in order to hold the section lift coefficient, c,, constant; 
tis the dimensionless span coordinate defined in Eq. (4); 
c, is the m.a.c.; and a and @ are the angle of attack and 
dihedral distributions over the wing in steady flight. 
Both a and @, as they occur in a,, @,, and a,, are evalu- 
ated for the deformed wing, so that ¢;,,°, ¢,,°, and ¢,,° 
do not correspond to a rigid wing, for the a and £ to be 
used will generally include the effects of “,, ao, and 4, 
both directly and indirectly (through elastic deflection). 

While the rotational motions p, g, and r produce 
small moment coefficient distributions, they may gen- 
erally be neglected for those wings where aeroelastic 
effects will be of importance. Hence, the section mo- 
ment coefficient may be written 


(36) 

The results [Eqs. (33)-(35)] are correct only for uni- 
form angular velocities and must be modified for har- 
monically varying velocities in accordance with the last 
paragraph of the Appendix. 


Cm = Cio + Cm 


METHODS OF SOLUTION 


There are several methods of attack which may be 
applied to effect a solution to the aeroelastic equations 
either in the form of Eqs. (11a) and (11b) or Eqs. (14a) 
and (14b). , 

The simplest possible configuration is that of a un- 
Swept wing (y = 0) of constant chord and stiffness, 
in which case Eqs. (lla) and (11b) are independent 
differential equations with constant coefficients and 
may be solved by the usual methods. 

If the chord and stiffness exhibit variations of the 
form (a + b£), (a, b, and c will be different for the 


AEROELASTIC 














SWEPT WING 183 





PROBLEM FOR 





chord and stiffness distributions) and the wing is un- 
swept, Eq. (11b) may be reduced to Bessel’s equation. 
This case has been treated by Hildebrand and Reissner.* 
The most important practical configuration is that of a 
linear (c = 1) chord distribution and quartic (c = 4) 
stiffness distribution, where the stiffnesses are propor- 
tional to the fourth power of the chord. In this partic- 
ular case, the transformation suggested in reference 2 
reduces to an Euler transformation,* wherein the in- 
dependent variable ¢ is introduced via the transforma- 
tion ¢= exp (1 — 7¢€), where (1 — 7¢€) represents the 
chord distribution. Fortunately, the last transforma- 
tion also reduces the simultaneous differential Eqs. 
(11) to differential equations with constant coefficients 
without imposing the restriction of zero sweep. Hence, 
it is possible to solve Eqs. (lla) and (11b) exactly fora 
large class of practically important problems. How- 
ever, in view of the many assumptions on which the 
solution rests, it will generally be more advisable to 
utilize approximate methods of solution because of their 
flexibility and simplicity, although an exact solution of 
Eqs. (lla) and (11b) might be useful in checking the 
approximate solution at one or more points. 

Equivalent forms of an energy solution may be ef- 
fected in several ways. The principle of least work* 
may be used, it being required that the difference be- 
tween the elastically stored energy (U) and the work 
(W) done by all the external forces remain stationary 
with respect to arbitrary first order variations about the 
true solution of the problem. The energy terms may 
be obtained from Eqs. (14a) and (14b), the unknown 
functions # and yg may be expanded in sets of deflection 
functions h, and ¢,, each satisfying the appropriate 
boundary conditions in Eqs. (13a) and (13b), and 
(U — W) may be made stationary with respect to each 
of the expansion coefficients to obtain a set of simul- 
taneous equations for these coefficients. 

In the case of an unswept wing, an alternative ap- 
proach to the simultaneous equations of the last para- 
graph is to formulate the value of g (dynamic pressure) 
for which the solution diverges. This value is evidently 
given by setting y = 0 in Eq. (14b), dropping the 
terms on the right side of the equation, and solving for 


q, Viz.: 
l —1 1 9 
0v\? 
[me fi cds | | f @( “) as | (37) 
0 0 Os 


It may be demonstrated that this critical value of g 
(proportional to the square of the divergence speed) is a 
minimum (not merely stationary) with respect to first 
order variations about the true solution for gy, this 
result being the Rayleigh principle’ for the second 
order, self-adjoint differential equation obtained by 
setting y = Oin Eq. (11b). By expanding ¢ in a set of 
functions ¢,, each of which satisfies Eq. (13b), and 
minimizing g, with respect to each of the expansion 
coefficients, a set of simultaneous equations for the 


Je | = 


| y¥=0 
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determination of these coefficients is obtained. The 
rapidity of convergence or relative accuracy of the 
solution may be determined by comparing the values of 
q. obtained by using an increasing number of terms in 
the expansion for ¢ or by using different approximations 
for the ¢,,. 

Perhaps the most direct method of approximating 
a solution to the differential Eqs. (11) is an extension of 
Galerkin’s method,® as exemplified by Bleakney’s ap- 
proach to the three-dimensional flutter problem.” The 
functions # and ¢ are expanded in sets of deflection 
functions hf, and ¢,, as in the previous approaches. 
These expansions are then substituted in the differen- 
tial Eqs. (11), Eq. (lla) is multiplied by /,,ds and Eq. 
(11b) by ¢mds, the resulting equations are integrated 
over the span, and a set (since m is arbitrary) of equa- 
tions analogous to Eqs. (12) isobtained. Since it is both 
expedient and advisable to require the expansion 
functions to satisfy the boundary conditions Eqs. (13), 
the stiffness terms may be integrated by parts to trans- 
form the set of equations to a set resembling Eqs. (14). 
Since m may take all values assumed by n, this set is 
sufficient to determine the expansion coefficients. The 
simplest application of this approach, wherein only one 
term is included in each expansion, is given in the next 
section. 

An iteration solution may be carried out by neglect- 
ing yg and y in the aerodynamic terms on the left side of 
Eqs. (11) and integrating the remaining stiffness terms 
to obtain the (static loading) approximations ¢° and 
y°. These approximations may then be substituted 
in the aerodynamic terms on the left, and the stiffness 
terms may again be integrated to obtain the next 
approximation. This procedure may be repeated until 
the agreement between two successive approximations 
is regarded as satisfactory. Such a method is recom- 
mended by the Army® for the calculation of reversal 
and divergence speeds and is also used in reference 2. 
However, in the case of a wing having appreciable 
sweepback, the process may not be convergent. In the 
case of the unswept wing (positive eccentricity), it 
may be demonstrated that the process converges for 
q < q., so that divergence of the process indicates that 
the wing divergence speed has been exceeded.* It has 
been shown further (by H. R. Lawrence) that a neces- 
sary and sufficient condition for convergence of the 
process, as applied to the more general case of the swept 
wing, is g < del, even though g, may be negative fora 
wing of appreciable sweepback and therefore have no 
physical significance. Now, in the case of no sweepback, 
q < q- is a necessary design condition, but in the case 
of a sweptback wing g may exceed ‘ae! if g, is nega- 
tive. Accordingly, the iteration method may gen- 
erally be expected to be unsatisfactory for swept 
wings. 





* By an appropriate normalization, the process may be made to 
converge at g = q, (see references 2 and 8). 
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INDUCTION EFFECTS 


In all of the approximate methods discussed in the 
last section, it is possible to include induction effects 
by replacing the products my and my with the actual 
(including induction effects) lift coefficient distributions 
corresponding to the geometrical angle of attack dis. 
tributions y and ¢, respectively. These lift coefficient 
distributions may be determined by the solution of the 
lifting surface integral equation, the method of Weiss- 
inger® appearing the most expedient for the swept wing, 
while the Prandtl lifting line solution’ will be sufi- 
ciently accurate for the case of the unswept wing. It is 
then necessary to solve the integral equation for each 
of the expansion functions. 

Hildebrand and Reissner? have calculated the sym- 
metrical divergence speeds of various unswept wings 
by an iteration process, coupled with a solution of the 
Prandtl lifting line integral equation. They conclude 
that divergence speeds predicted for typical wings on 
the basis of strip theory, using the two-dimensional lift 
curve slope (m), are from 17 to 40 per cent too low. 
However, if their results are further examined, it ap- 
pears that the error is considerably reduced (to less than 
5 per cent in all practical cases) if the three-dimensional 
lift curve slope, O0C;,/Oa, is used in conjunction with 
strip theory. The latter lift curve slope was calculated 
for the rigid wing in the comparisons made by the 
writer, but in practice it might be obtained from wind- 
tunnel measurements on the wing alone. 

Further investigation of the unswept wing shows that 
the asymmetrical divergence speed may also be pre- 
dicted with satisfactory accuracy (5 per cent or better) 
by calculating m for a wing of half the true aspect ratio. 
In fact, it appears that, with the suggested evaluations 
of m for symmetrical and asymmetrical distributions, 
all aeroelastic properties that do not involve the control 
surfaces may be adequately predicted by the use of 
strip theory. On the other hand, control surface char- 
acteristics of an unswept, elastic wing could not be 
satisfactorily predicted on the basis of strip theory 
alone (compared to lifting line calculations). However, 
in view of the success of strip theory in predicting the 
divergence speeds and other aeroelastic characteristics, 
it appears that its failure in predicting aeroelastic char- 
acteristics may be due primarily to the incorrect dis- 
tribution of ¢;;° specified by Eqs. (27) and (29). If so, 
Cis° may be correctly evaluated, but strip theory may 
be used in calculating the additional lift distribution 
due to elastic deflection. 

In the case of an unswept wing, ¢,,° may probably be 
calculated with sufficient accuracy on the basis of lifting 
line theory (as it was in the aeroelastic results cited 
above), although calculations should be made on the 
basis of lifting surface theory to check this prediction. 
Unfortunately, neither the Prandtl lifting line theory 
nor the Weissinger Method? will suffice to compute ¢," 
for a swept wing, and lifting surface theories such as 
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FORMULATION OF 


those of Falkner''! (for subsonic speeds) or Lagerstrom!” 
(for supersonic speeds) should be used. 

Whether or not strip theory, together with an appro- 
priate value of m (determinable on some simple, a priori 
basis), will be adequate for predicting the aeroelastic 
characteristics (other than control surface) of a swept 
wing remains to be seen, but it appears from the results 
for the unswept wing that such is probably the case. 


EXPLicitT SOLUTION 


An explicit solution may be obtained from the energy 
Eqs. (14) by writing 


g(t) = Ag(&) see y (38a) 
h(é) = BIf(£) ese v (38b) 


[the sec and ese y are introduced for later conveni- 
ence in writing a(£), cf. Eq. (1)] where f(£) and g(é) are 
assumed to satisfy Eqs. (13a) and (13b), respectively. 
In principle, Eqs. (88) may be regarded as exact, but in 
practice f and g will generally be estimated, and the 
resulting solution will be approximate. Introducing 
Eqs. (38) in Eqs. (14), together with the dimensionless 
span coordinate & = (s//), yields 


1 f' .. (dg\? ‘genes 
GJ \|—] d& — qml cos* ec*g*dé | X 
1 Jo of 0 
rd 
A sec y + | om sin y cos? f ec* -. g ae| xX 
0 f 


ad 
Bese y = qi cos® vf c(ec)’ + Cm)gdéi + 


if ok? x ( (“eae if of (=) eat (39a) 
-| amr cos? vf cfg at sec +(} f EI X 
(=!) dé + qmil* sin y cos ¥ . gta |B Bese y = 
o£ 0 ar 
gl? cos vf ccf de +22 foe (ys f dt (39b) 
where g in (2/g) 


g) and (4/g) is the acceleration of gravity 
and not the torsional deflection function. 

Eqs. (39) may be directly solved for A and B, thus 
completing the solution for ¢ and h. In order to dis- 
cuss further the various aspects of the aeroelastic 
problem, it is, however, convenient to introduce a 
number of dimensionless parameters. (In order to 
simplify the discussion, the inertia loading terms will 
be omitted in the following, although they could be in- 
cluded in a similar manner.) Thus, Eq. (39) may be 
cast in the convenient, nondimensional form 


(\—-1A+pB=C, (40a) 
—A+(k+1)B=C, (40b) 
A = (Gu/Cun)(4GSo/meogS* cos vy) (41) 
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Ex SEI, 
«= (=)\(— eae (42) 
Cx/\m A.R. sec? y gS? sin 


a= (Ci2/ Cir) (43) 
jie (Cy ‘Cx) (44) 


1 
Gu = (S/cob)? f GJ (Og O£)*dé (45) 
0 
eo 
Ex = (S/cob) f EI (0°f/d€?)2dé (46) 
1 
Ci td Cul, = aym) f ec?[c;,° + (Cm/e) |g dé (47) 
0 


1 
C2 = CnC. = (1/m) [ éc,°f dé (48) 
0 


1 
Cu = f ecg? dt (49) 
0 


1 
Cy = f ec*g(Of/OE)dE (50) 
0 


1 
Cy = f fg dé (51) 
0 


1 
Co f Gf(Of/d#)dt (52) 
0 


(18) = (18/110) 
[Jo = £100) | 


where the bar indicates the value of a quantity relative 
to its root section value* and the subscript zero indi- 
cates the root section value. The solution of Eq. (40) 
is given by 


(53) 


= [« + 1)G — 2G] A" (54a) 
= [A — 1) + »JAa- (54b) 
A=(A-— Ll(k +1) + ww (55) 


DIVERGENCE SPEED 


The critical value g, (or divergence speed) is obtained 
by equating the determinant (55) to zero, which yields 


h = 1 — w(l + «)7 (56) 


Since \, and «x, are both inversely proportional to q,, 
Eq. (56) does not represent an explicit solution for q, 
except in the special case x, = © (corresponding to 
infinite bending stiffness or zero sweep). 

Consider now the case of a wing with positive eccen- 
tricity (e.a. aft of l.a.c.) and positive sweepback, so 
that both \ and « are positive for positive values of q. 
Remembering that torsion of such a wing produces 
washin, whereas bending produces washout, it is phys- 
ically obvious that divergence is impossible for suffi- 
ciently low values of EJ and therefore of x, and it 
follows that \, < 0, or ww > 1. Conversely, if it is 


- © Except i in the case of C; and ( 








486 JOURNAL OF THE AERONAU 
assumed that the wing is sweptforward (eccentricity 
still positive), it is physically evident that a positive 
value of u, must exist for sufficiently small values of x, 
and it follows that uy < 1. Hence, if it ts assumed that 
the deflection functions are identical for positive or nega- 
tive sweep, it follows that 


pv = | 


and Eq. (56) may be rewritten 
— (*) at -(2\(2 ‘\G R. sin *) (GJ)o 
aad Ke - Cu Fx» 2€ cos* Y (El )o 


(58) 
It is evident that even moderate sweepback of a wing 
with appreciable aspect ratio may make divergence im- 
possible, whereas sweepforward may produce a low 
divergence speed. 


AERODYNAMIC DERIVATIVES 
An approximation suggested by Eq. (57) is 
Fa) 
Ci _ uC» = 7 IC, 


It is found to be accurate when c,, 
Cm is appreciable, as in the case of control surfaces, but 


(59) 


= (in C,, while if 


constant over the range of integration of Ci, Eq. (59) 
may be written 

C, = uC, = rv IC, (60) 

Tr = 1 + (€p/ec,’) (61) 


Substituting Eqs. (57) and (60) in Eq. (54), it is found 
that A and B are suitably approximated by 


A = C2(q) = eul'C.22(q) (63a) 

= €,0(q) = (eu T)—'C,9,(q) (63b) 

2 = xA-! = [A — (1 — e})]-' = & = (64a) 
1 


( 
% = NA! = [x + (1 — &)] 7! = € 'Q, = (64b) 
(<) = (\(2 ‘\( 2eo cos* Y \(=*) ” 
e={- > (69) 
PN Gu C. A.R. sin ve GJo 
The angle of attack and dihedral distributions are ob- 


tained by substituting Eqs. (38) and (63) in Eqs. (1) 
and (2), now including a and 8°. The results are 


a(t) = a%(&) + [g(&) — (eT) oo (66a) 
= w(t) + [(eu)g(t) — f’'() CM (66b) 
B(E) = BE) + [tan yg(é) + (eu ’-! cot or (&) )1C,24(q) 
(67a) 

= 6) + [(exI’) tan y g(t) + cot ¥ f’(E)]E.22(q) 
(67b) 


For large values of euI’ it is evidently more convenient 
to use Eqs. (66a) and (67a), whereas for small values of 
this parameter Eqs. (66b) and (67b) are appropriate. 

The lift and moment coefficients may be calculated 
by substituting Eqs. (38) and (63) in Eqs. (18)-(25). 
The results may be written in the form 


(C,/Cxr°) = 1 + Rp,Q(q) 
= 1 — k,9(q) 


(68a) 


(68b) 
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(Cap/Cy®) = 1 + Ry Qi(Q (69a) 
= 1 — ky&(q) (69b) 

(C,/C,) = 1 + R,Q(q) (70a) 

= 1 — &,2(q) (70b) 

ky, = (mC,/C1,)[Czr, — (eel) —C,,] (71a) 
kr, = (mC2/C,°)[C.. — (euxT)Cz] (71b) 


—mC i -F La... 
Rwy, = ( ager )4( Ji L, — (ex) “tC, + 
‘3s Vu Cm 


b tan +\,... 
( tc uv, — (exI) ‘Cunlt (72a) 


<Cm 


—mC2\ { (a + Xm 
oa, = ( = y( ic. — (eul’)Cz,] + 
C Mu Cm 


(° = Vic  ~ Pte ; (72b) 

i - + (mc u, — (wT)—Cy,] (73a) 
i-s (Jc ue — (@T)Cu, (730) 
Cr, = (e0b/: sf fg dé (74a) 

Cr, = (c0b/: sf c(Of/O&)dE (74b) 

Cu, = (eb S) f exe (75a 

Cu, = (d/S fe R(of/deedE — (75b 


DEFLECTION FUNCTIONS 


The accuracy of the foregoing results depends on the 
basic assumptions and on the deflection functions f and 
g which are selected. The simplest possible nontrivial 
functions satisfying Eq. (13) are the deflection functions 
for a uniformly loaded cantilever beam and are given by 
2§ — # g(1) = 1 (76a) 


g(é) = 


f(&) = (1/3) (GE — 48 + &), = f(1) = 1 (76d) 


Some improvement in accuracy is obtained by using the 
static deflection functions of the wing. Assuming the 
chord and stiffness distributions 

(1 — 7é¢), €l)hjewi—=]|-T=r (77) 


é(é) = 


GJ(é) = El(é) = &(&) = (1 — Te)‘ (78) 
the deflection functions obtained by substituting 
OF /Os ~ é(€) andO7/O0s ~ c?(é) in Eq. (5) and integrat 
ing are 


37°g(é) = In &(&) + 3)[E(—)-3 — 1] (79a) 





~_ 


WET. papa 


(69a) 
(69b) 
(70a) 
(70b) 
(71a) 


(7 1b) 


72b) 


73a) 


3b) 


4a) 


Ab) 
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6T*f(E) = (1 + 37? — 7%) + c() In &(é) — 
., +1... fA... 
i+5T 3% r? 1 e(&) — 9 é(é)—! + 
73 
(5)ew 


The set (79) has not been normalized to unity at the 
tip, as in the case (76), but otherwise Eqs. (79) reduce 


(79b) 


to Eqs. (76) for T = 0. 


EXAMPLE 


In order to illustrate the results of the foregoing 
sections and to bear out the important consideration 
that aeroelastic effects must be given in the design of a 
swept wing, a numerical example will be used to con- 
clude the paper. 

Consider a tailless aircraft 


characteristics: 


having the following 


A.R. sec? y = 12 

7, = c(l)/c(0O) = 0.5 
GJ(é) = EI(&) = &(&)3 
MG marS?/ (GI )o = MG mareS*/(ET)9 = 10 (maximum) 
e= 0.15 

50 per cent span, 20 per cent chord control surfaces 


The panel aspect ratio is chosen as large as is generally 
feasible, and the ratio mqS*/(G/J)o = 10 corresponds 
roughly to what might be expected for a rather heavily 
loaded wing designed for adequate gust factors. 

The various dimensionless coefficients will be calcu 
lated for the simple deflection functions (76).* (Fur- 
ther calculations have shown that the error in using 
Eqs. (76), compared with Eqs. (79), is about 5 per cent 
in the important parameters uw and x.) Substituting in 
Eqs. (45)-(52) and Eqs. (74) and (75), for é(&) = 1 


(£2), the results are: 


Gy = 0.373 Gn/Cu = 1.79 
Ex = 1.750 Fx2/C2 = 5.58 
Cy = 0.208 wm = (Ci2/Ci,) = 1.59 
Ce = 0.330 vy = (Cn/Cn) = 0.672 


II 


Cy = 0.211 
Cy = 0.314 
Ci, = 0C)/Oay = 0.298 
C,, = OC2/0ao = 0.256 
Ci, = 0C,/06 = 0.175(—0a/085)qI'5 cos ¥ 
C,, = 0C,/06 = 0.203(—0a 06), COS ¥Y 


Ci, = 0C,/0(pb/2U) = 0.165 
C,, = O0C2/0(pb/2U) = 0.176 
Ci, = 0.568 
Cy = 0.934 


Cn, => 0.340 
Cw = 0.519 
C,, 1.43 


19 
* Actually, g(t) = ('/2)(3& — &%) was used in these calculations 
in order to compare with certain other results; however, the use 
of g(t) = 2 — £2 would not appreciably affect the results. 
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Cra, = 0.814 

Cis = 0.842(—da/05)_I"; cos y 
O, = 0.646(—0a/05),, cos ¥ 
Ci, = 0.794 

C», = 0.560 


2p 
It is seen that uv = 1.07, checking reasonably well with 
Eq. (57). Using the experimental results of 7.R. 721," 
it is found from Eq. (61) that Ty; = —1.4, it being as- 
sumed that both c,,, and c°), exhibit the same variation 
with sweep angle and Mach number, although there is 
some considerable doubt about the last assumption. 

Substituting the above data in Eqs. (41), (42), and 
(65), it is found that 


r == 1S sec ¥ q . Gg = (q q mar.) 
0.372 ese y G7 
0.078 cot ¥ 


K 


lI 


€ 


Three different configurations, obtained by giving 
the above wing different angles of sweep, will be con- 
sidered. The results shown in Table 1 are then self- 
explanatory. 

The aeroelastic corrections to the various aero- 
dynamic derivatives may now be obtained by sub- 
stituting the foregoing data in Eqs. (68)-(70). The 
results for the pitching moment derivatives will be 


evaluated for (a/c»,) = 0.25 (see Table 2). 
The reversal speeds are given by 
G\ Cita = 9 25(>¢-) 0.764 0.780 
gi Crs = 0 1.83 0.757 0.730 
g\ Cus = 0 —7.40 3.34 
gC = 9 2.06 1.75 2.08 


It may be observed that the C,, also exhibits a finite 
reversal speed for the swept wings. It would appear 
that this phenomenon is probably not physically 
possible and is due to the error in the assumed deflection 
functions; however, if this indicated reversal speed is 
not too closely approached, the resulting error in C,, 
is probably small. The aeroelastic reduction in rate of 
roll is given by 


p/p? = (C/ Cy") (Cup"/ Cy) (50) 
and, from the foregoing data, is given (approximately) 
by 

I I] III 
p/P’ 1-049¢g 1-049¢9 1-0.38¢ 


for speeds not too close to reversal. 

Comparing the results for the swept and unswept 
wings, it is evident that the most serious aeroelastic 
problems introduced by sweepback are loss of static 
stability (reduced C,,,) and trim effectiveness (reduced 
C,, and Cy,). The latter may be particularly serious, 
since the increased control surface deflections may 
mean an increased drag. The problem of aileron rolling 
effectiveness is, on the other hand, no more serious in 
the presence of sweep, although this result may be by 
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: TABLE 1 ' 
Quantity Eq. Ref I II III a 
¥ 0 45° 63 
N (41) 1.77 q°! 6.75 97! 10.597? 
x (42) x 0.525 q7! 0.417 gq"! 
€ (65 . 0.078 0.040 
eu : 0.124 0.064 
euls ; — 2 —0.174 —0.089 
A = (1 — e~!) (58) l —11.8 —24.2 
qe A = 2X. 4.4 —0 573 —( 435 
FLY (71a) 0.813 ; roe 
2a (71b) 0.702 0.730 
EM ax (72a) 0.813 = 
Mea (72b) Bio 0.702{1 + 0.322(cm/a)] 0.730[1 + 0.261(cm/a)) 
C1x3° (20) 0.416 15, 0.416 Cis° 0.416 Cl5q 
Fig (71a) ~1.61 . 
126 (71b) 1.61 1.53 
Cas° (21) —[(a/cm) + 0.835] Cxrg° —[(a/cm) + 0.688] Cz5' 
(72b) d (a/c 3.8: l¢ 
ku, 72b 0.620 Bes 3.10(a ; m) 0.581 Sex 3 83(a/cm) 
1+ 1.20(a/cm) 1 + 1.45(a/cm) 
Ci5° (22) —0.153ci5° —0.153¢15° —0. 153159 
Fig (73a) —1,31 en 
125 (73b) ae 1.22 1.16 
Ci, (22) —0.139m —0.139 m —0.139 m 
an, (73a) 0.970 a cas 
i, (73b) 0.960 1.00 
Q; (64a) [4.77 q-! — 1]7! o pace 
f (64b) (0.525 g-! + 0.922)—} (0.417 g~! + 0.960)! 
TABLE 2 
I II IT! 
Clas —— “ 
—, 1 + 0.170 g(1 — 0.21 9) 1 — 1.34 9(1 + 1.75 @)~ 1 — 1.75 q(1 + 2.30 9)— 
Clay 
Cm a a . — a . ; : 
= s 1 + 0.170 g(1 — 0.21 g)=! | — 3.07 g(1 + 1.75 9) 1 — 3.58 g(1 + 2.30 9)! 
Mag 
(2 - 0.25) 
Cm 
Cig : , ss » as 
Gai 1 — 0.338 g(1 — 0.21 g)7 1 — 3.07 g(l + 1.75 9) 1 — 3.67 g(1 + 2.3 9)! 
Lb 
Cg : = ; 
Cus 1 — 1.615 g(1 + 1.75 g)= 1 — 2.00 G1 + 2.3 g)- 
(4 = 0 25) 
Cm 
Cg ” = ~ 
Gai 1 — 0.274 g(1 — 0.21 9)" | — 2.32 g(1 + 1.759)" 1 — 2.78911 + 2.39)7! 
t 
Ci 
Ga 1 + 0.203 g(1 — 0.21 g)7' | — 1.83 gl + 1.759) 1 — 2.40 g(1 + 2.3 g)= 
1 
Pp 
no means universal. In addition, the aeroelastic upwash at the section; and mp, the section lift curve 


modification of the dynamic stability derivatives for a 
swept wing may be expected to be different, and prob- 
ably more important, than for an unswept wing. 


Appendix 


Lirt DISTRIBUTION ON A RIGID SWEPT WING 


Consider a chord section, c’, perpendicular to the 
l.a.c., and let U’’ be the relative wind velocity compo- 
nent along this section when it is at zero incidence; U’, 
the component (directed outboard) along the (swept) 
l.a.c.; W, the upwash at zero incidence; w, the actual 


slope, possibly corrected for compressibility and induc- 
tion effects but not sweepback. The section lift is then 


] - 
= c’ ( pu”) m( =) 


If the section is rotated about the swept l.a.c. to an 
incidence a’ and about its neutral axis to a dihedral 
(with respect to the swept axis) 8’, the resultant upwash 


given by 
dL 


7 ds 


l (Al 


is given by 
w= U'a’ — U"p’+W (A2) 


The section lift due to these rotations is then given by 


m/a)} 


Q)~! 


Li 


an 
ral 
sh 


bo 
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Eq. (Al) as 8’ = —asin y + Bcos y (A4e) 
pL’ nia wae : Substituting Eqs. (A4) in Eq. (A3) yield: 
1 = c! — m(U'a’ — UB’ + W) (AB) intial niin taal 
lL = c’(pU?/2)(mpo cos ya (A5) 
In the case of steady flight with forward velocity U, whence the effective section lift curve slope is given by 
incidence a, and geometric dihedral 8 (all angles as 
m = My COS ¥ (A6) 


sumed small), 

If the result [Eq. (AG)] is taken for the sweepback cor- 
rection and mp is presumed to be the section lift curve 
U sin 7 (A4b) slope, uncorrected for induction effects and compressi- 
bility effects, a corrected three-dimensional lift curve 
slope at subsonic speeds may be obtained from the 
j acos y + 6 sin y (A4d) formula (empirical and speculative) 


U’ = Ucos y (Ada) 
8 had 
W=0 (A4c) 


~ 
ll 


__ Mo COS y(1 — M? cos? y)~' 


m= — - 
Mo COS Y " z eee ae (A7)* 
1+ (™ ) (1 — AF cos* ») "(1 — M*)~"? (1 + 7) 
wv A.R. 
where the factor (1 — MW? cos? y)~'* is a compressi- / = c’(p/2)mo[(U + u — ry) cos y + 
bility correction on 79, being the usual Prandtl-Glauert'® (v— ré) sin y] [((U+u—ry)at+ 
correction based on the component of flow normal to the e ; 
Sanco ‘ , ; v — rz) w ny y, 
la.c.; (1 — M?*)~ “isa correction to give the effective \ B+w+ py+qi] (All) 


aspect ratio; and 7 is a correction to account for the Dropping products among (u, v, w, p, g, 7) as second 
deviation of induction effects from those produced by _ order yields 


an elliptic lift distribution.’ ss ead . 
fe . L=c'(p/2)l OS 2u — 2ry 
The result [Eq. (A7)], for 4 = 7 = O, has been cle mola cos y(U + 2n yr 
checked by Weissinger® against the results calculated vtan y — ré tan y) + Bcos y(v — rz) + 
the basis of his lifting line theory and was within 2 
on the basis of his lifting li eory ithi con vie + on + el (ALD) 


per cent for a wide range of sweepback angles and taper 


ratios. Weissinger also cites supporting evidence for The section lift coefficient may then be evaluated as 


the use of half the actual aspect ratio in Eq. (A7) in the ! 
case of asymmetrical lift distribution. 6 a (mo cos ¥) ‘ck + 

The available 7 factors!” have, of course, been calcu- c’(pU*/2) 
lated on the basis of lifting line theory for unswept 2u v 

ag a, : . -~-+ —tan y — — (tan y + 2y)]}+ 
wings, but their use for moderately swept plan forms is i I ’ 
probably sufficiently accurate to warrant inclusion. : : 

' , ; Sa fv rt w py. qt\ ; 

For large sweepback, + must be obtained empirically. a( 7; —7 ) + 7 7 + rf (A13) 


In the case of small disturbed motions of the airplane 
with a given C.g. position, the resulting values of U"and if mm, cos ¥ is replaced by from Eq. (A6), differentia- 
U" at the point y on the l.a.c. are given by the expres-  ti9n yields 


sions 
Cy, = OC,/O0(u/U) = 2ma (A14) 


U’ = (U+u— ry) cos y + (wv — rz) sin y (A9a) 
Cy = 0¢,/O0(v/U) = m(atan y + 8) (AI15) 


U" = (U+u — ry) sin y — (v — r&) cos y (A9b) 


_ j Cy = Oc; O(w U) = m (A16) 
W =wtpyt+ gt (A9c) 
The results [Eqs. (A14) to (A16)] are due to the veloc- 
tC = Xm ) ti (A10) . . 
: oF ta TTY ity components on the l.a.c. In addition, the rotational 


where (u, v, w) are the velocity increments and (p, g,7) motions p, g, r produce nonuniform chordwise velocity 
are the angular velocities about the (x, y, z) axes, respec- distributions. These distributions may be treated by 
tively. It is assumed that the c.g. is essentially on the thin airfoil theory in the same manner as camber. 
Substituting Eqs. (A4d), (Ade), and (A9) in Thus, if (c/2) cos @ represents the chordwise coordinate 
and w(6@) represents the vertical velocity distribution 


(wing relative to fluid), the section lift and moment 
(with respect to the quarter chord) may be written" 


x axis. 
Eq, (A3) yields 


* Note: Eq. (A8) deleted. 
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C= -2 f [w(0)/U](1 + cos @)d6 = (Al7) 


Cm = off [w(6)/U](cos 20 — 1)dé6 — 5 (A18) 


For a rotation of angular velocity w about the quarter- 


chord axis 


w(0) = —(we/2)[cos 6 + (1/2)] (A19) 

whence 
C, = 2r(wce/2U) (A20) 
Cm = —(2r/4)(we/2U) (A21) 


The results [Eqs. (A20) and (A21)] have been obtained 
without reference to sweepback and for an implicit lift 
Hence, they should be multiplied 
The net 


curve slope of 27. 
by (m/2m), and c must be replaced by c’ 
angular velocity of a chord section about the l.a.c. is 


given by 


w = —psiny + qceos y + r(asin y — Bcos y) (A22) 
whence Eqs. (A20) and (A21) yield 
mc Cos” ¥ ; 
C= 50 [—ptany +q+r(atan y — B)] (A23) 
mc COS” ¥ 
Cm = [ptan y — gq — r(atan y — B)] (A24) 


16U 


Taking the derivatives of Eq. (A23) and adding them 
to the corresponding derivatives of Eq. (A13) yields 


Oc; E Ros | \25 
C, = -~ =m — —- Sn + cos AZ) 
ly d(pb 2U) $ b Y Y ( 


OC;  ( # ) a ee 2 A926 
+ = = mM “0s? (A26) 
_ O(gem/ U) Cm 2 Cm a ee 
dc; | =) i + p+ 
, = 2s ~— an y ) 
4 O(rb/2U) — = © 


B 
2Ea + (‘) cos? y(@ tan y — a)| (A27) 
)? 


The section moment coefficient are obtained by differ- 
entiating Eq. (A24) with the results 


OC m m ‘) _ 49 

Cmp = - = sin 2 (AZS) 
: 0(pb/2U’) 16 \6 LU 

OC » m ( “) F 19 

Cmg = a Sa =) oo8" 4 (A29 
¢  O(gem/U) 16 \Cn 

OC m m (<) : 
Cm, = = — cos* y(@ tan y — 
d(rb/2U) = \,) cos’ ve Y B) 

(A30) 


The effective angle of attack distributions listed in 
Eqs. (30)-(35) are then obtained by dividing Eqs. 
(A14)-(A16) and (A25)-(A27) by m. The moment 
coefficients given by Eqs. (A28)—(A30) can generally 
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be neglected for wings of reasonable panel aspect ratios 
(A.R. sec? y). 

The foregoing results are based on constant perturba- 
tion velocity components (u, v, w, p, g, 7). In stability 
calculations, however, perturbations exhibiting a har- 
monic time dependence are of more general interest, 
and it is necessary to take nonstationary aerodynamic 
effects into account, even though the frequency be small. 
These effects have been considered elsewhere,'> and 
if f is the frequency, Eqs. (A20) and (A21) must be 
modified to read 


we 2U 
C¢, = 2x ( ) [ 2077216 — In ( : )| (A351 
21 mh 


—(mr/2)(we/2U) A32 


Cn = 


Accordingly, for harmonic motion, the terms propor 
tional to cin Eqs. (A25)—(A27) must be multiplied by the 
term in square brackets in Eq. (A31), while Eqs. (A2S 

(A30) must be doubled. It should also be added that 
the compressibility correction for stability derivatives 
is not simply the Prandtl-Glauert factor (1 — ./*) 
but is considerably more complex, even for linearized 


“/3 


theory and low frequencies. 
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Landing-Gear Oscillations Due to Unstable 
Skidding Friction 


J. E. WIGNOT* anno FREDERIC M. HOBLITt 
Lockheed Aircraft Corporation 


ABSTRACT 


It is shown theoretically that unstable (self-excited) skidding 
oscillations of an airplane landing gear may occur under service 
conditions. The primary requirement for such oscillations to 
occur is that the coefficient of sliding friction between the tire and 
the runway decrease with increasing skidding velocity; such a 
variation ordinarily exists for a rubber tire on wet pavement 
Consideration is also given to the influence of other variables 
on both the initial stability of the oscillations and the maximum 
amplitudes, loads, and stresses that may be reached. 


INTRODUCTION 


— RECENT YEARS there have been a number of 
heretofore unexplainable landing-gear service 
failures that have occurred while the airplane was skid- 
ding on a wet runway. Reference | indicated that a 
possible cause of such failures was the occurrence of 
self-induced oscillations of the gear. Since this phe- 
nomenon may occur at relatively low speed, it was easy 
to overlook it as the real cause of the failure and at- 
tribute the failure instead to an accumulation of the 
effects of previous hard landings. The purpose of this 
paper is to present a quantitative theory to explain 
such self-induced oscillations. 

The landing gear of an airplane must be capable of 
fulfilling two functions: first, absorbing the energy 
associated with the vertical velocity of descent of the 
airplane and, second, absorbing the energy associated 
with the forward velocity of the airplane. The first is 
accomplished with the aid of a “‘shock strut,’’ which is 
usually a hydraulic type damper, while the second is 
accomplished by braking the wheels and converting the 
into an elementary  solid-friction type 
“damper.’’ The discussion that follows will be con- 
fined to the latter function. 

Inasmuch as the characteristics and effectiveness of 
the braking system depend directly upon the coefficient 
of friction between the tire and the runway surface and 


inasmuch as investigations of coefficients of friction of 
3 


airplane 


tires on various surfaces under various conditions* 
indicate that the coefficient of friction varies with the 
velocity of slip of the tire, it is of interest to investigate 
the effects of such a variation on the drag forces ex- 
perienced by the landing gear. 


Presented at the Session on Structures—Dynamic Loads, 
Seventeenth Annual Meeting, I.A.S., New York, January 
24-27, 1949. 

* Structures Engineer 

t Stress Engineer. 
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THE COEFFICIENT OF FRICTION 


The factors that appear to affect the coefficient of 
friction are many and varied. Some of the more ob- 
vious ones are runway surface material and finish, 
runway condition (wet, dry, oil film, etc.), tire tread, 
tire inflation pressure, direction of skid, and velocity of 
skid. The actual variation of the coefficient of fric- 
tion with these factors is complex and little understood, 
but all investigators agree that in many cases, partic- 
ularly on wet surfaces, the coefficient of friction de- 
creases with increasing velocity of skid; such a varia- 
tion in coefficient of friction will be assumed for this 


investigation. 


FUNDAMENTAL PROBLEM 


For simplicity, the following possible service condi- 
tion will be considered. An airplane is moving straight 
forward on a level runway with constant vertical load 
on the gear and with the brakes locked so that the 
wheels are skidding. Some disturbance is then en- 
countered, such as a rough spot in the runway, which 
imparts motion to the landing gear. (The problem will 
remain essentially unchanged if the disturbance is 
simply the original locking of the brakes.) The prob- 
lem is to determine the nature and magnitude of this 
motion and the magnitude of the drag loads which the 
landing gear experiences as a result. 


ELEMENTARY THEORY 


The motion of the landing gear will be studied with 
the airplane considered fixed in space while the runway 
appears to move aft with a velocity Vo. The compli 
cated physical system may be idealized as shown in 
Fig. 1. 

The effective mass .\/ of the wheel and lower strut 
is assumed to be concentrated at the axle and is at 
tached to the lower end of a flexible strut that at its 
upper end is rigidly attached to the airplane. 

The following horizontal forces act on the mass V7: 

(1) A spring force provided by the flexible strut, 
applied at the axle and equal to Ax, where A is the 
spring constant of the strut and x is the horizontal dis- 
placement of the axle. 

(2) A friction force uW, applied at the ground, where 
u is the coefficient of friction. Since the coefficient of 
friction is a function of velocity, it must be expressed 
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Fic. 2. Assumed variation of coefficient of friction. 


in terms of the variable velocity of the landing gear, 
dx /dt. For small variations in velocity, the coefficient 
of friction is given approximately by 


uw = wo t+ S {(dx/dt)[(h + r)/h]} 


where wo is the coefficient of friction at velocity Vo, S 
is the slope (positive as shown in Fig. 2) of the curve of 
u versus velocity, and the factor (i + r) ‘his introduced 
to account for the fact that, with motion about an ef- 
fective center of rotation A (Fig. 1), the velocity at the 
ground differs from that at the axle. The friction force 
is then 

h +7 wo 


W = oWV Ss 
“i r + h dt 


(3) A force due to structural or other damping, 
assumed to act at the axle and to be proportional to ve- 


locity. It is given by c(dx/dt), where c is the damping 
constant. 

(4) An “inertia force,’ ./ d*x,;dt*, acting at the 
axle. 


The other forces acting are the reacting force and 
moment at the top of the strut and the upward ground 
force W. 
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Making use of d’Alembert’s principle (whereby an 
acceleration may be replaced by its equivalent “‘inertia 
force’ and the laws of statics applied to the resulting 
force system) and applying a virtual displacement 
consisting of a bending of the strut in which the locked 
wheel moves about an effective center of rotation A 
through an angle 66, the following work equation re- 
sults: 


Act+r 


= _adx 
Kx(h60) — E +3 W - 


|e + r)d0 + 
h dt 


dx d*x 
Cc F (h60) + MV iP (héd) = 0 


at at 


Rearranging, dividing by /60, and letting S[(h + r)/h}? 
= S’ gives 


7, d*x 1 SW dx K W h+r 
d ie = ) 2 xX = po ( 
dt? dt _ h 
The general solution of this linear differential equation 
is: 

x =e ™(C, sin Bt + C. cos Bt) + x, (2) 
where 


a = (c — S'W)/2M 
8 = V(K/M) — [(c — S’W)/2M]? 
X, = (uoW/K)[(h + r)/h] 


It can be seen that x, is the static displacement of the 
axle in the x direction under the action of a constant 
frictional force woW acting at the point of contact with 
the runway. 

The constants C, and C, are obtained from the initial 
conditions. If it is assumed that the plane is skidding 
with the brakes locked and without oscillation of the 
landing gear and if a small displacement A from the 
occurs on account of a 


Il 


equilibrium position x = Xx, 
bump, rough spot, or similar irregularity in the runway, 
then att = 0,x = x, + Aanddx/dt = 0. From these 
two initial conditions, the constants are found to be 


C,= A 


Cy = (a/B)A 


Eq. (2) now becomes, if these values are substituted 
and the equation put into nondimensional form by di- 
viding by x,, : 


oS 
=l1+-—e — sin 6t + cos Bt (3) 
< Xe B 


If the initial disturbance were to be provided by sud- 
denly applying the brakes while the airplane is coasting 
freely, the initial displacement would be given by 


A= —xX;, 


and Eq. (3) would be modified to 


x : 
-=l]—e* E sin Bt + cos a | (4) 
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LANDING-GEAR 


Examination of the terms on the right-hand side of 
Eq. (3) leads to the following conclusions: 

(1) The factor [(a@/8) sin Bf + cos Bt] will be a 
steady periodic function with a radian frequency equal 
to 8and a maximum amplitude equal to V1 + (a/B)?. 

(2) As indicated by the first term, any oscillation 
that does occur will occur about the static displacement 
as an equilibrium position. 

(3) The factor (A/x,Je~™ modifies the oscillatory 


function 


[(@ 8) sin Bf + cos pt] 
in such a manner that the two curves given by 
1 = (A/x,)V 1 + (a@/B)2e7% 


provide an envelope within which this oscillatory func- 
tion must be contained (Fig. 3). Thus, according as 
the value of a is positive, zero, or negative, the ampli- 
tude of the oscillation will decrease, remain constant, 
or increase, respectively. 

(4) Thus, the stability of the oscillations can be 
determined by examining the value of a. 

(5) Since a is given by (c — S’W)/2M and since WV 
must be a positive quantity ‘c and W will likewise be 
positive, as will S if the coefficient of friction decreases 
with increasing velocity as assumed), it follows that 
the stability criteria are: 

(a) Ifc — S’Wis positive (damping predominat- 
ing, or S negative), then the oscillations are 
damped. 

(b) Ifc — S’W = 0, the oscillations are stable. 

(c) Ifc — S’W is negative (runway friction vari- 
ation term S’Il’ predominating), then the 
oscillations are unstable. 

(6) The amplitude of the oscillations at any instant 
is directly proportional to the magnitude of the initial 
disturbance, although it will be independent of the 
source of disturbance. (Note, however, that with 
unstable oscillations, if a certain amplitude is reached 
after a certain number of cycles following a given initial 
disturbance, then with a much smaller initial dis- 
turbance the same amplitude will be reached even 
though only after a greater number of cycles.) 
EXTENSION OF ELEMENTARY THEORY TO NONLINEAR 
EFFECTS AND MAxXImuM LOADS 


An incorrect conclusion that might be drawn from the 
elementary theory is that, once an unstable situation is 
encountered, the amplitude of oscillation would actu- 
ally increase indefinitely with time. 

Such a conclusion is indicated by the elementary 
theory because of the assumption of a linear variation 
of both damping force and friction force with velocity. 
While such an assumption is an excellent approximation 
for small oscillations, as the amplitude of oscillation in- 
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Fic. 4. Nonlinear variation of damping and friction forces. 


creases it is necessary to consider the effects of the non- 
linear nature of these force variations. 

The nonlinear nature of these forces is illustrated by 
Fig. 4, in which the effective force F at the axle is 
plotted as a function of skidding velocity V between 
the tire and the runway. 

Curve A in Fig. 4 is the curve of friction force, u.W 
[(2 + r)/h], in which the factor (4 + r)/h corrects the 
force uW at the ground to an effective force at the 
axle. (The effective force at the axle is related to the 
actual force at the ground in such a manner that the 
two forces do an equal amount of work in a virtual ro- 
tation about the center of rotation A, Fig. 1.) 

Curve B in Fig. + gives the damping force, which is as- 
sumed to act at the axle. At a skidding velocity equal 
to the airplane velocity |p there is no relative velocity 
between the landing gear and the airplane, so the 
damping force is zero at the point V = Vo. At smaller 
or larger values of V, dx/dt assumes positive or nega- 
tive values, respectively, and the damping force varies 
accordingly. The displaced coordinate V’, against 
which the damping force is plotted, is not equal to 
dx/dt, however, but instead is given by (dx/dt)[(h + 
r)/h\. 

Curve C is the algebraic sum of curves A and B. 
This curve is the key to the important characteristics 
of the nonlinear motion. 
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NONLINEAR STABILITY CRITERIA 
Since for infinitesimal values of dx/dt or V’ (Fig. 4) 
curves A, B, and C are all approximately straight lines, 
the elementary theory applies, and the stability criteria 
given for the linear system are correct also for the non- 
linear system. 


MAXIMUM AMPLITUDES FOR UNSTABLE OSCILLATIONS 


When an unstable mode of vibration has some initial 
excitation, the amplitude of vibration increases with 
time as long as the net work done on the system is posi- 
tive. When the net work done on the system is nega- 
tive, energy will be removed from this mode of vibration 
and the amplitude of vibration will decrease. There- 
fore, the criterion for the maximum amplitude of vibra- 
tion is that the net work done during a complete cycle 
shall be zero: 


¢¢42e = 0 


where E is work done on the system and § indicates 
the integral over one complete cycle. 

Returning to Fig. 4, a horizontal line F = F% is 
drawn through curve C at | = Vp. Forces measured 
from this line (to curve C) will be denoted by F’. Then 
the total work during one cycle is given by 


R= § Fdx 
= ££ (Fy + F’)dx 
= ft Fo dx + f F'(dx dt)dt 


But § Fy dx is the work done throughout a complete 
cycle by a constant force and is therefore equal to zero. 
Also, dx ‘dt = [h/(h+7r)|V’. Consequently 


E = [h/(ht+n|g F’V' dt (5) 
dE = [h/(h + r)|F’'V’ dt 


II 


if only that portion of dE is considered which does not 
integrate to zero over the cycle. 

Since dt is always positive, dE is positive when F’ and 
V’ are both positive or both negative, and dE is nega- 
tive when F’ and V’ are of opposite sign. 

Accordingly, at small values of V’ as indicated by 
amplitude (a) in Fig. 4, dE is positive throughout the 
cycle, since curve C within this range is at all points 
situated so that V’ and F’ are of the same sign. At 
larger values of |’ (both positive and negative) as occur 
during parts of the cycle of amplitude (b), F’ changes 
sign while V’ remains the same sign, so that dE becomes 
negative. As the amplitude increases, more and more 
negative work will be done during the cycle, until at 
some velocity amplitude the net work will be zero. 
This will be the maximum velocity amplitude. 

To determine the maximum amplitude associated 
with any given curve C, a graphical or numerical pro- 
cedure can readily be devised, based upon Eq. (5), for 
evaluating the total work per cycle, /, for any given 
velocity amplitude. (It may be assumed in the evalua- 
tion of F that the frequency will be equal to the un- 
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damped natural frequency, w = VK/M and that dis- 
placement and velocity both vary sinusoidally 90 
out of phase with each other; such an assumption js 
probably accurate at maximum amplitude, which is the 
condition being investigated.) E will be evaluated for 
various velocity amplitudes until one is found for which 
E = 0. This is the maximum velocity amplitude, 
Vmar., from which the maximum displacement am- 


plitude of the axle is obtained by the relation 


X maz, = (Vmaz.’w)[h/(hk + r)] (6 


where X = x — Xs. 

It should be emphasized that, for a system in which 
unstable oscillations occur, the oscillations will ap- 
proach asymptotically the same limiting amplitude, 
no matter how small the initial disturbance may have 
been. However, a certain period of time is required 
for the oscillations to build up before they closely ap- 
proach this limiting amplitude; consequently, it may 
sometimes happen that before the limit is closely ap- 
proached a change occurs in external conditions, such 
as runway surface or airplane speed, which would make 
the system stable and thus halt the increase of ampli- 
tude. For example, since the airplane is skidding with 
brakes locked, its speed might decrease to a region of 
smaller S or to zero before the oscillations had time 
to build up. If such a phenomenon occurs, it is appar- 
ent that the smaller the initial disturbance, the smaller 
would be the maximum amplitude reached in the time 
available before the conditions change. In practical 
design, this phenomenon probably need not be consid- 
ered, since the critical skidding condition would occur 
with the large initial disturbance due to sudden ap- 
plication of the brakes. 

It should also be noted that the limiting velocity am- 
plitude depends only upon the friction and damping 
characteristics of the gear and not upon its stiffness or 


mass. 


MAXIMUM OSCILLATORY LOADS AND STRESSES 


The maximum load induced by the oscillations is 


: : ; ee & 
| here a X mar K = J mar -_ 
h+rao 
h h 7 
1 ae : M , = ] Fai ‘ egy 4 7 } 
“A+r a hay VKM \ 


Thus the load is decreased by a decrease of either A or 
M, the other being held constant. It may be either 
decreased or increased by a decrease in w, depending 
upon whether / or KA is held constant. Furthermore, 
the load is directly proportional to Vimaz.’; 
many cases it will be approximately proportional to 
” reaches Vo, dE 


hence, in 


the airplane speed, Vo, since, when | 
suddenly becomes negative and extremely large. 
Changes in M and Vyaz.’ will have the same effect 


on stress as on load. When K is changed, however, 
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LANDING-GEAR OSCILLATIONS 


the strength changes as well as the load, so that the 
effect on stress is not immediately apparent. If the 
“spring” is stressed in bending (as in Fig. 1) and A is 
varied by varying wall thicknesses, with depth of sec- 
tion and all other geometrical features held unchanged, 
then stress in the ‘‘spring’’ is proportional to displace- 
ment, as given by Eq. (6). (If the “‘spring’’ were to 
carry load by means of axial rather than bending stress 
and A were varied by varying the cross-sectional area 
of the axially loaded members, stress in the “‘spring”’ 
would still be proportional to displacement, and the 
following results would still apply.) For these condi- 
tions, it follows from Eq. (6) that 


w © Veer’ V M/K (8) 


. , , 
i mar a J mar. 


Thus the stress in the “‘spring’’ is decreased by an in- 
crease in w, Which may be brought about by decreasing 
M or increasing K. The same qualitative conclusion 
would be reached if, in varying A, geometric similarity 
of the cross section instead of constant depth of section 
were to be maintained, although Eq. (8) would then 
not be exactly correct. 

Parts of the landing gear which would not be altered 
in the process of varying A would encounter a stress 
proportional to load, as given by Eq. (7), so that, for 
these parts, stress as well as load would be decreased by 
a decrease of A. Such parts would probably include, 
for example, wheel bearings, drag strut, and other at- 
taching parts. 


OBSERVATIONS OF UNSTABLE SKIDDING OSCILLATIONS 


The unstable skidding oscillations described by the 
foregoing theory have been observed on an instru- 
mented airplane skidding on a wet runway, in tests 
described in reference 1. For this particular airplane, 
installation of an hydraulic drag strut damper, which 
had been designed primarily for the purpose of reducing 
dynamic overloads due to wheel spin-up, was found to 
reduce greatly the maximum amplitude of the skidding 
oscillations. 

Unstable skidding oscillations that have been ob- 
served qualitatively in other engineering structures are 


described in reference 4. 
CONCLUSIONS 


(1) Landing gears may experience large oscillatory 
drag loads under service conditions due to the physical 
phenomenon of a coefficient of friction that may de- 
crease with increasing velocity of slip. 

(2) The magnitude of such drag loads is independ- 
ent of the source of the initial disturbance which starts 
the oscillation and independent of the magnitude of 
the disturbance provided all conditions, including air- 
plane speed, remain constant for a sufficient period of 
time. 

(3) An increase in the equivalent static load 
W[(h + r)/h]? on the gear, or in the slope of the fric- 
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tion coefficient vs. velocity curve (positive as shown in 
Fig. 2), increases the probability of the oscillation oc- 
curring, increases the rate at which it builds up, and 
increases (perhaps only slightly) the maximum, oscil- 
latory drag load. 

(4) An increase in the damping present decreases 
the probability of the oscillation occurring, decreases 
the rate at which it builds up, and decreases (perhaps 
only slightly) the maximum value of the oscillatory 
drag load. 

(5) An increase in mass (stiffness remaining con- 
stant) increases the drag loads and stresses throughout 
the landing-gear structure. 

(6) An increase in stiffness (mass remaining con- 
stant) Such 
stiffness ordinarily involves an increase in strength of 
the parts being stiffened, so that the drag stress in these 
The drag stress in all other parts 


increases the drag loads. increase in 


parts is decreased. 
of the landing gear is increased. 

(7) An increase in frequency of oscillation (as 
produced by variation in stiffness or mass or both) 
decreases the displacement amplitude of oscillation, 
may increase or decrease the drag loads, and ordinarily 
decreases the drag stress in parts that affect the stiff- 
ness. 

(8) Before the theory can be applied quantitatively 
to predict susceptibility of a particular design to un- 
stable skidding oscillations, it will be necessary to in- 
vestigate further the variation of coefficients of friction 
with velocity of slip. Such investigation would be 
particularly desirable at high speeds, in order to deter- 
mine whether the trend toward higher landing speeds 
makes the landing gear more or less susceptible to failure 
due to these oscillations. 

(9) Before the theory can be applied quantitatively 
in the early stages of design it will be necessary to find 
a means of estimating with sufficient accuracy the mag- 
nitude of the structural and other damping. (After an 
airplane is constructed, the damping constant can 
readily be measured. ) 

(10) 
stable skidding oscillations is to introduce large damping 
into the system in the form of an hydraulic damper. 


A possible practical means of control of un- 
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Relations Between the Characteristics of 
a Wing and Its Reverse in Supersonic Flow 


A. H. FLAX* 
Cornell Aeronautical Laboratory, Inc. 


SUMMAR\ 


On the basis of linearized supersonic aerodynamic theory, cer- 
tain relations between the aerodynamic characteristics of a given 
wing and another wing having the same plan form but with lead 
ing and trailing edges interchanged are derived. The latter 
wing is called the reverse of the first. Only the case of ‘‘purely”’ 
supersonic flow (Mach wave always behind wing leading edge 
is considered. It is shown that in such cases the slope of the 
lift curve, the drag due to lift, the wave drag at zero lift, the 
damping in roll, and the damping in pitch are the same for any 
given wing and its reverse. In addition, a simple relation be 
tween the static pitching moment derivative of a wing and the 
lift derivative due to pitching of its reverse is obtained 


INTRODUCTION 


— PROBLEM OF CALCULATING the lift and drag of 
supersonic wings according to the linearized theory 
has been dealt with by a number of authors.'~‘ In 
general, they have been concerned with the solution of 
the problem for wings of some -particular plan form. 
In this paper, some relations between the aerodynamic 
characteristics of a wing and another wing having the 
same plan form but with leading and trailing edges 
interchanged are derived. The latter wing is called the 
reverse of the first. Only the case of ‘pure’ super- 
sonic flow (Mach wave behind the wing leading edge) 
will be considered. In particular, it is shown that, if 
the wing is inverted by interchanging the positions of 
the leading and trailing edges, the lift due to angle of 
attack, drag due to lift, and drag at zero lift are un- 
changed. In this way, each solution already obtained 
gives the solution for a second case in addition. Further, 
it is shown that the damping in roll and pitch are 
governed by similar relations. Finally, a relation be- 
tween the static pitching moment derivative of a wing 
and the lift derivative due to pitching of its reverse is 
obtained. As one immediate application of the results, 
the case of the arrow or ‘‘delta’”’ wing is easily treated 
by consideration of the reversed arrow wing on which 
the lifting pressure at every point is the two-dimensional 


value. 
Basic EQUATIONS 


The linearized equation for isentropic, nonviscous 
supers nic flow is, using Puckett’s notation,! 

Received July 18, 1948. Revised and received March 30, 
1949. 

* Head, Aerodynamics Department. 
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ee» beet 82 ee l) 
Ox" oy? ds? 

where 

M = U/a 

u = 0¢/Ox uo = O+4u 

v = 0¢/dy vy =v 

w = 0¢/0zs w= w 


Here u’, v’, and w’ are the fluid velocities in the direc- 
tions x, y, and z, respectively; U is the free-stream 
velocity directed along the x axis; and a is the velocity 
of sound. The velocities u, v, and w are assumed to 
be small compared to U, thus representing small per- 
turbations in the main free-stream flow. The wing 
surface is assumed to lie in the x-y plane, and the wing 
is assumed to be extremely thin. Puckett! has shown 
that the solution for the potential ¢ may be written as 


If w(t, n)dé dn 
(x,y) = — a = ——— : z) 
rV (x — £)? — By — n)? 


where 8 = VM? —1 , and where, for the case where the 
Mach line lies behind the wing leading edge, the inte- 
gration is carried out over all points, &, 7, on the wing 
lying in the fore Mach cone of the point x, y, as shown 
in Fig. 1.1 In effect, the solution, Eq. (2), represents 


+ Note that this equation differs in sign from that given by 
Puckett. This merely represents a change in the direction of 
positive w and is accounted for in calculating w from the boundary 


conditions. 
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A WING AND ITS 
the potential produced by a distribution of supersonic 
sources at the points &, n, of strength w (&, ). 

The increment of pressure in excess of free stream 
pressure produced at any point is 

P=-pUu=-jU0¢, (3) 

to the accuracy of the linearized equations, which 
leads to the pressure coefficient 


Cp = P/('/2)p U2? = —2¢6,/U (4) 


The integral of Eq. (2) is finite in spite of the fact 
that the integrand becomes infinite at the edges of the 
Mach cone. If, however, an attempt is made to calcu- 
late , by differentiation under the sign, it is found that 
difficulties arise, since the resulting integral is infinite. 
Hadamard? has dealt with just this problem in his work 
on the general solution of the Cauchy problem for hy- 
perbolic partial differential equations. He introduced 
the concept of the finite part of an infinite integral of 


the type encountered. Thus, according to Hadamard, 
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| fs gw, 0) dn dé 

o(x, vy) = f f win at (2a) 
Sad V(x — 8)? — By — 0)? 

where the symbol denotes the finite part. This 


integral may be differentiated under the sign, and the 
resulting integral will still be defined by its finite 
part. In performing the differentiation with respect 
to x, attention must be given to the variable upper 
limit. This may be dealt with in the usual manner, 
leading to 


oe w(x, y) _ 
r\¥s J — 


B 
If w(t, n) (x — §) dn dé (5) 
» 
m[(x — £)? — By — »)*) 


This expression will be used in the present work. 

The method of evaluation of the finite part of an 
infinite integral is given in the Appendix. An ap- 
proach to the problem not employing Hadamard’s 
method is also included in the Appendix. 


RELATION BETWEEN LIFTS 


If A is the area of the wing, the wing lift coefficient due to angle of attack may be written 


Cr ss 


} ¥ Ww t ’ if ~ a ¢ ™ " 4 wx, Vv) 
SS SS (w(t, n)/U)(s £) dé dn fe dy 4 SS (v, eae (6) 
TA A [(x — £)? — B?(y —_ n)?] 8 A A B U 


where the integration with respect to x and y ranges over the entire wing. 
The reciprocal relation between a wing and its reverse may now be derived. 


For ease of visualization, this is 


illustrated in Fig. 2 by an arrow wing and a reversed arrow wing. The derivation is by no means limited to this 


case. 
from the boundaries of the wing plan form as &, 7. 


on the wing, and let the fore Mach cone trace of the point x’, y’ bound an area C’ on the reverse wing. 
seen that these are physically the same segment of area relative to the boundaries of the plan form. 


Consider some point, £, 7, and a point x’, y’ on the reversed wing which is located at the same distance 
Let the aft Mach cone trace of the point £, » bound an area C 


It will be 
The 


increment of lift produced on C by the source at point &, 7 is then 


iC, = 
der - 


4dé dn if | (wi, n) /U) (x — 8) dx dy + 
cJ [(x — §)? — By — n)*]"” 


w(t, n 


) 
48.0 dé dy (7) 


Returning to the reversed wing, consider the increment of lift coefficient produced at point x’, y’ by sources at 


points £’, »’ lying in the area C’. This is 


dC,’ = 


For both wings at a constant angle of attack a, w(x, y) 


is a constant 


w(x, y) = Ua 
or 
w/U=ea (9) 


An examination of Eqs. (7) and (8) for dC; and dC,’ 
discloses that, since the integrals are symmetrical in 
x,y, and &, », except for sign which is compensated by 
the change from fore to after cone, dC, = dC,’. In 
other words, the increment of lift produced on the aft 


ABU 


tdx’ dy’ w(t’n’)/U)(x’ — &’) dt’ dn’ . 4w(x’y’) . 
a i . 
7A cS [(x’ — &')? — By’ — 9’)?]” 


Mach cone area C by a source at &, 7 is equal to the in- 
crement of lift produced at the corresponding point 
x’, y’ on the reversed wing by a distribution of sources 
of the same strength per unit area distributed over the 
area C’ in the fore Mach cone. 

The total wing lift coefficient as given by Eq. (6) 
is obtained by integrating Eqs. (7) and (8), with re- 
spect to dé dn and dx’ dy’, respectively, over the entire 
area of the wing. Since there is a one-to-one corre- 
spondence between é, 7 and x’, y’, it is readily seen that 
C, = C,’. Further, since the coefficient of drag due 
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FIG. 3 


to lift is simply equal to C,a, it is seen that the drags due 
to lift for the two wings are similarly equal. 

This result may be applied directly to the arrow 
wing. Puckett! has found by direct integration that, 
for an arrow wing with the Mach wave swept behind 
the leading edge, the lift coefficient is given by 

Cr = fa/V M? — 1 
which is exactly that for a two-dimensional unswept 
supersonic airfoil; it was pointed out that this is sur- 


It was, in fact, this result that led to the 
It will now be shown that this 


prising. 
present investigation. 
result may easily be obtained from the results of the 
study thus far. 

The lift for a reversed arrow plan form is readily ob- 
tained by a consideration of Fig. 5. If the Mach 
waves lie behind the trailing edge, as is always the case 
for ‘‘pure’’ supersonic flow on the arrow counterpart, 
then the fore Mach cone intercepts an area on the wing 
which is always the same as if the point were on an in- 


finite two-dimensional wing. This is easily seen by a 

















(7) 


FIG. # 


comparison of Figs. 3a and 3b. Thus, the pressure at 
any point, which is a function only of the source dis- 
tribution in the fore cone, will be the same as if that 
point were on an infinite two-dimensional wing. There- 


fore, for the reversed arrow, 
Cr = 4a/V M? — 1 


Consequently, by the reciprocal relation, this is also 
the lift on the arrow wing. 

Fig. 4 illustrates some typical pairs of wings to which 
the reciprocal relations may be applied. 


RELATION BETWEEN WAVE DRAGS AT ZERO LIFT 


Again using Eq. (5), the expression for the wave drag coefficient at zero lift of a supersonic wing may be written 


: } “ff f fF rx) A(E)(x — £)dE dn dx dy 4 
Cp = : — ar 
TA A [((x — §)? — By — »)7]” A 


dx dy (10 


ji 
A B 


where } is the local slope of the airfoil measured in the xz plane and the range of integration with respect to x and 


y is the entire area A. 


Returning to Fig. 2, the increment of drag coefficient produced on the area C by the source at point &, 7 is 


4dé d 
dCp = = = 


TA 


[(x — &)? 


SS. 


— p(y — n)*] 


A(x) A(E)(x — £) dx dy 4 


+ MEd Ed (11 
5 AB 7 


For the reversed wing, the increment of drag coefficient produced at point x’, y’ by sources at points £’, n’ lying it 


the area C’ is 
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é 4dx' ae A(x f= 2 ) dt’ d 4 
dCy’ = = Ew aot 2— 8% ee — aa + AB d?(x") dx’ dy’ (12) 


Again it is found that the integrals given by Eqs. (11) and (12) are symmetrical in x,y and &, 7 except for sign, 
and again they are equal for the same reasons given for the lift integrals. Thus, the increment of drag at zero 
lift produced on the aft Mach cone area C by the source at &, 7 is equal to the increment of drag at zero lift pro- 
duced at the corresponding point x’y’ on the reversed wing by the distribution of sources over the area C’ in the 








fore Mach cone. 

The total drag coefficients for the wing and the reversed wing are also equal by the same reasoning as was used 
in the proof of the equality of the lifts. This proof of equality of wave drag at zero lift applies somewhat more 
generally than the proof of equality of lift, in that it holds even where the Mach wave lies ahead of the leading edge, 
if the top and bottom surfaces are symmetrical, so that there is no flow from top to bottom. 

One immediate consequence of these results is to make all the data for the zero lift wave drag of ‘‘delta’’ wings, 
given by Puckett,’ valid for the corresponding reversed arrow wings. This is of considerable value, since the zero 
lift wave drag of a reversed arrow wing is by no means so simple to calculate as its lift. 








nr 
» 
‘ RELATIONS FOR DAMPING IN ROLL 
The damping moment, /, due to a rolling velocity p, may also be calculated by the use of Eq. (5). In this case 
w(t, n) = pn (13) 
The damping moment derivative C;, is defined by 
7 
Cy = 1/[(C/2)p U?Ab(pb/2U)] (14) 
where bis the span. The value of the a derivative in roll is then given by 
pb. ny(x - =— - td E di n dx dy . 
*Cip PF ial dy (15) 
20°” A a5 | (x — t)? — 6? (y — n)2]°? + 
In Fig. 2, the increment of rolling moment a on the area C by the source at point £, 7 is then 
pb _ Ap dé dy ha ny( x - ¢) dx dy 4pn? 
sure at dC, = + —— didn 16) 
mage 0° ~ De ot c[(x — 8? — By — m3)" | ABUb \ 
if that | On the reversed wing, the increment of rolling moment produced at point x’, y’ by the sources at points £’, n’ lying 
There- | inthe area C’ is 
a 4p dx’ -| (x’ — g’ dé’ d n’ . 4py”? " 
7 dC,! = co + PX dx’ dy’ (17) 
at TA Te” — £')? — B%(y’ — 9)?]'* | ABOb 
is also | 
whic By reasoning similar to that used previously, dQ, = 7” y : 
F dC,’ and C, = Cy’. — §) d§ dy 
This result may readily be applied to obtain the ‘\ cla — 5? By ; yay7* + 
damping in roll of an arrow wing. First, the damping : M 
in roll of a reversed arrow wing is established. For the PO cll (x — £) dé dn 
reversed arrow wing, the pressure coefficient at any dill 2py 90 
vritten point when the wing has a rolling velocity p is (referring Ly c[(x — £)? — B%Xy — »)?]'” + UB (<U) 
to Fig. 5) ; , , , 
, The second integral in this expression may be shown to 
(A } 2pm be zero by geometrical reasoning. To each point &, » 
(x — £) dé dn opy on the right of the line 1-1 in Fig. 5 corresponds a point 
) x and “\ | f- — = $ 7 : equally distant from the line on the left, having the 
[(x Fe 21 B same value of (x — £) and a value of (» — y) which is 
” (18) equal but opposite in sign, thus canceling the effect of 
Now ; the first point so that the total value of the integral is 
(11 ‘ow may be written as zero. The sum of the first integral and the constant 
. n=y+(9— 9) (19) term in Eq. (20) is nothing more than the pressure at 
ying in ’ a point on a two-dimensional wing having an angle of 
and attack py/l’—which is, of course, the local angle of at- 
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/ which, when integrated becomes 
7 4/4 4 5 l = TT? b° Co fon 
} L = oS had ee (22) 
Ble, VJ in — 2 4 12 
/ 
vy where Cy is the wing chord at the center. This leads to 
> 
C, = 1/3V M2 - 1 (23) 
This value then is also valid for the arrow or delta wing 
by the relations obtained previously. 
X 's RELATIONS FOR DAMPING IN PITCH 
. 
The damping moment due to a uniform rotational 
FIG. S& velocity, g, about some line in the plane of the wing, 
perpendicular to the direction of flight, say x = 0, 
tack produced by the roll. Thus the local pressure Will be studied next. The local velocity at any point 
at every point on a rolling reversed arrow wing is the 1S 
same as that on a two-dimensional wing having the w(t, n) = gt (24) 
same local angle of attack as that produced by the roll- 
ing velocity. The pitch damping derivative may be defined as 
The damping moment in roll for the reversed arrow M C 
wing may then be written as Cnq = (25) 
ee ee ee * (/) p02 A C (qC/0) 


where C is the mean chord of the wing. Applying Eq. 


(5), the derivative for damping in pitch is 


ydA (21) 


gee. p (; . 0") Yi 
pV — 1 02" -0/2 
, Ex(x — £) dé dn dx dy 4 , 
Cn = —S a — x? dx dy 
Fal c [(x — §)? — BxXy — )2]"” ~ BAC? 


The increment of pitching moment derivative produced on area C by the source at point £y is 


( { t(x—® dxdy 4 
Sf. a 5 a = E° dé dy (27 
| [(x — 8)? — BxXy — n)2]”" ~ BAC? 


On the reversed wing, the rotation is assumed to occur, and the moment is calculated about the sameaxis. This 
means that on the reversed wing the axis of rotation is the trailing edge, say XY’ = Cy. The increment of pitching 
y’ by sources at points é’, n’ lying in the area C’ is 


(26) 


—4d id 
dCnq = 
wAC? 


moment derivative produced at point x’, 
fas (Cy = x’)(Cy = #/)(x" — #/)d8" dn’ 
c [(x’ — £’)? — By’ — 9’)?]'” 


about the leading edge due to rotation about the leading 


ctr 


re es tdx’ dy’ 
Cum = wAC? 


— —— x'*dx'dy’ (28 
BA C? 


There are the following relations between the coordi- 


nates on the two wings: edge is 


(29) 4 Pe ie / x? 
M=--— pU? [q=dA 
VM? -—12 ( 


Cg’ by the same line of 
The necessity 


, ; 
Co —i = é, Co = é = xX 


Thus, dCng = dCng’ and Cy. = 
reasoning as has been used previously. 
for choosing the same axis on the wing and its reverse 
should be noted, however. , 
The damping moment in pitch for an arrow wing will M = -— D 
next be calculated by using the reverse arrow wing, for 
which the calculation is much simpler. Returning to 
Fig. 3, it may be seen that the pressure at every point 
on a reverse arrow wing in pitch will be the same as 


which becomes, when integrated, 


This leads to 


Cue = -—-S 3V M2 —_ 1 (31 





that on a two-dimensional wing having the same local 


Thus, the moment for the damping derivative of the reverse arrow wing | 


z velocity component due to pitch. 


(22) 
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A WING AND 
for rotation about its leading edge and for the arrow 
wing for rotation about its trailing edge. 


RELATION BETWEEN THE STATIC PITCHING MOMENT 
DERIVATIVE AND DERIVATIVE OF LIFT DUE TO PITCHING 


The static moment coefficient for a wing may be 
related to the coefficient of lift due to pitching for the 
reversed wing ina simple manner. The lift derivative 
due to pitching about a given axis may be defined as 


-f fe 4ddédn 
A TAC 
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Crug = L/{(*/2) p 0? A(QC/0)] (32) 


The static pitching moment derivative may similarly 
be defined as 


Cua = M/[(*/2) p U2 CAa] (33) 
Referring again to Fig. 2, the static moment derivative 
on the wing on the left measured at x = 0 (a line 


through the apex) is 


acl S, tdé dy (34) 


For the reversed wing on the right, the lift derivative in pitch, for pitching about the same axis (the trailing edge) 


is 


SS, 4dx’ dy’ 
A 7wAC 


where Cy is the maximum chord of the wing. If it is 
noted that Cy — &’ is equal in value to x, then, since the 
remainder of the inner integrand in both expressions is 
symmetrical in x,y and &,n, 


CL,’ = Cita and (36) 


’ , 
Cig = Ca 


According to these results, the static derivative of 
moment of a wing about any axis is equal in magnitude 
to the derivative of lift due to pitch about the same axis 
of its reverse wing. It should be noted that the axis is 
to be held fixed in the wing 1ather than relative to the 
flow. 

Eq. (36) may be applied to find both C,, and Cyq 
for the arrow wing, since the calculations for the re- 
versed arrow wing are again particularly simple. 


CONCLUSION 


For wings in supersonic flow, when the Mach wave is 
swept behind the leading edge, it has been shown that 
for sections symmetrical about the chord the lift and 
the drag at zero lift remain the same if the positions of 
the leading and trailing edges are interchanged. In 
the case of the drag at zero lift, this conclusion holds 
even when the Mach wave lies ahead of the leading edge 
as long as the upper and lower surfaces are symmetrical. 
The damping derivatives in roll and pitch, relative to 
axes fixed in the body, have also been shown to be 
equal for the wing and its reverse. Further, it has been 
shown that the static moment derivative of a given 
wing is equal to the lift derivative due to pitch of its 
reverse wing, when the moment and the pitching mo- 
tion are taken about the same axis, fixed in the wing. 
It should be kept in mind that these conclusions are 
subject to the usual limitations of the linearized theory 
of supersonic flow. 


fee. (Co — es 





x’ — ) d&! dr 
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nn staf ff ste ay 


The concept of the effect of a source on the area in its 
after cone, introduced in this paper, may be useful in 
the solution of various supersonic flow problems, 
where it may lead to simpler integrals for certain cases. 
In addition, it appears that the methods of this paper 
may be applied to establish relationships between a 
wing and its reverse for several stability derivatives 
not considered in this work. 


After the present work had been completed, it was 
called to the author’s attention that W. D. Hayes had 
arrived at similar results. These results are given 
only for lift and drag and were obtained by different 
considerations from those used in this paper. 
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Appendix 


The loading integrals utilized in this paper, which are infinite and require the use of Hadamard’s finite part ofan | 
infinite integral, may be avoided by the use of a less direct approach to the problem. This will be demonstrated, 
and it will be shown that the same results are obtained with Hadamard’s finite part of the loading integral. 

Consider first the increment of load produced on the point x,y by a uniform distribution of sources on a strip at 
a distance £ from the leading edge of the supersonic wing in Fig. 6a. This strip represents a fairly general case, 
since one end of the strip is on the edge of the fore cone of point x,y, while the other is on the wing boundary, which 
may be an arbitrary function f(£). The point on the fore cone leads to the singularity in the lift integral. The 
increment of lift produced at point x,y, may, however, be calculated without introducing this difficulty; this, in 
fact, is what has been done in most of the previous work on the linearized theory of supersonic wings. First the 
potential at point (x,y) due to the strip under consideration may be written as 


f(§) ee en 
ie = cn f dé dn/ V (x — £)? — By — 0)? 
7~(e-O/ 


for unit source strength. This is readily integrated to give 
f(é) 


d@ = (dt/rB)\are sin [B(y — n)/(x — &)] 
9 = te ~ 2/03 
or 


d@ = (dt/x@)(are sin { BLy — f(é)|/(w — &)f — (#/2)) 


The velocity due to this potential is 





dt Bly — fOVW/@ — 8 


do, = = = = 
mV 1 — Bt Ly — f(é)]/(e — &)}? 


The increment of load on an area dx dy at point xy is then 


— bly ~—fOV/e— §? | 
V1 — Bh [Ly — f(6)/(x — 8}? 


dP = — 2 dx dy dé 
To evaluate the total load on area dx dy due to the sources lying in its fore cone, it is only necessary to integrate 
with respect to &. This introduces no new difficulties, the singularity as & — x having been dealt with in the 
earlier part of the paper. 

Next the inverse case is considered. The problem of calculating an increment of load produced by a source at 
point £, » on the corresponding strip in the after cone of the reversed wing (Fig. 6b) requires a more careful treat- 
ment. Consider first the effect of a finite source strip at £ of coordinates » + m, and» — m. The potential at any 
point of the aft strip due to the forward source strip is 


"2 +m cant 
do = (dé/m) | dn/V (x — &)? — By — n)? 
”—m™ 
or 
7+rM™ 


d@ = (dt/r) are sin B[(y — n)/(x — &)] 


7; = & 
The increment of velocity is 
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The increment of lift produced by this velocity is 


ly - + m)] Bly — (n — m)) 
Ne dy 


a -_ af Cr __@&- st - 7 i? eF __@-8! 
‘ #2) ;= (n- 7 ay f(x) yi _ 62 k ~ Gn - xy) 


es it 


x—€& 
Integrating 
pUdé dx | ‘ , ‘| -€5 =< 
P= ————j| (x — §)* — ly — ) 
d (x — E)xB? (x — & 8? Ly (n + m)] e 


7 — §)? — #*[y — (2 — »)]' ‘] Po ¢ 7) 


I(x) 


a ap pU dé dx miKG » ¢)? - (t— = = nf] - |e =o t)? — oll —(n + ™)| | ; — 
(x — t) r 
on \2 _ &- 8) is . 2 2 - ™ 
oo = Be +m + | (x — &)? — 67 f(x) — (9 — m) ( 


It may be seen that if m becomes zero, the increment of load becomes zero, which would be expected. To obtain 
the effect of a small strip of width dy (n, = dn/2), it is necessary to take the derivative of dP with respect to n; at 
m = Oand to multiply this by dy/2. 


| _ &— ) _ 
; pUdédndx | B n| B° Lf(x) — (7 + ™m) )] 
dP = tim (x — Bas? SS 7 . “2 — - rer 
Ryne 2 Vic - 9 - af - (x D+ a (x — £)? — Bf (x ” m1? 
B 
x— §) 
| — + 

. ‘aw n| ; a*[f(x) — (a — mw) 


—@=8 an V(x — 2 — BLf@) — (n — mE 


(x — &)? — =| - 
\ 3 
The first and third term in the brackets cancel as 7; — 0, giving 


dP — p = dg dn dx — a. x) a ” — —<— 


ap = PU ded dy B(f(x) =») Vi - fees) 


rB (x — &)? x—€t 


This is seen to be equivalent to the value for dP developed for the direct problem. 
It may be of interest now to illustrate how Hadamard’s finite part of the integral leads to the same result. 


Consider the direct case: 
ait? f(E) 
pl (x — &) ) dn 
dP = dx dy dé f aa 
T y ~ (« - o/8 [(x — 8)? — BUy — »)2]'” 


This may be written 


Le eam — &) 
dP = _BU(x — £)dx dy dé I 8 J . 
7 f(e) [(x — &) — Bly — »)] “°[@ — &) + BY — 9)] 


Hadamard has shown! that the finite part of an integral of the form 


f A(x)dx/(x, — x)'” 
a 


which is appropriate to this problem may be written as 
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A(8n) = 1/[(x — §) + By — »)]' 


where, for this case 


and 
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Then 
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Further 
| 7.( b)'/*1( t) + 8 |’ | 
y—m— max -— £)"[@ —é B(y — m)]”’ 
= — f(®)1 | it m1 (2)" 2 m 
T= - — + lim - aT 
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The second term goes to zero in the limit. 


pU dx dy dé Bly — f(é)] f(é) 7? 
iP = — 
' 7B (x — £)? [\i- LE <3 


which is the same result as was obtained before. 


Using this result 
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Spanwise Lift Distribution for Sweptback Wings 


Alan Pope and William R. Haney, Jr. 
Georgia Institute of Technology, Atlanta, Ga. 
April 15, 1949 


HE METHOD OF Schrenk! for obtaining the spanwise load 
eee of an unswept wing has been discussed and 
extended by Flatt,? who (a) advanced a procedure for accounting 
for twist and flaps and (b) noted that the Schrenk-Flatt method 
required approximately 5 per cent of the time needed by other 
methods. From a limited number of calculations, the authors of 
this letter concur with Flatt. 

However, Schrenk and Flatt did not consider sweptback wings. 
The methods of Weissinger® and Theilheimer‘ are frequently used 
for this problem but are somewhat lengthy. The following 
formula, which assumes that the effect of sweepback on the non- 
dimensional span loading curve is linear, appears to give results 
closely approximating experiment while ré€quiring only a minimum 
of time, a many-point span loading curve being obtainable in 
about 1!/2 hours. 

The method is as follows: 

(1) Find the nondimensional spanwise lift distribution for 
zero sweep (cc;/@CL), = 9 according to Flatt-Schrenk. (c = local 
chord, c; = local lift coefficient, € = average chord, and C; = 
wing lift coefficient, A = angle of sweepback at quarter chord.) 
For an untwisted wing of constant airfoil section, (cc,/E€Cz), = 9 is 
simply the mean between the geometric chord and a half-ellipse 
having the same area and the same major axis as the wing. 


Reduce it locally according to the relation 


CCI cci \ 2y 
CCL A CCL ee b 


where y = local span station, ft.; and b = wing span, ft. 

Eq. (1) is demonstrated in Fig. 1, where it is seen that the sec- 
ond term represents a triangular loading diminution that is 
largest inboard. Subtraction of this loading from the unswept load- 
ing yields a new wing Cr which may be found by integrating the 


curve. 


[2(1 — cos A)] (1) 


Thus 
Cra = So cei/ECindn (2) 


where » = nondimensional span station 2y/b. 

(3) To obtain values for other lift coefficients, increase or de- 
crease the unswept nondimensional loading by multiplying it by 
the desired amount. Then correct as before to obtain the load- 
ing with sweep. 

Four comparisons of this method with others and experiments 
are given in the drawings. The following comments are perti- 
nent: 

Fig. 2. A.R. = 4.66, 31° sweep, taper ratio 2.26 to 1. The 
method of Eq. (1) closely follows both experiment and Weissinger. 
At the tip it is slightly closer to experiment than Weissinger. 
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Fic. 2. Spanwise distribution of section lift coefficient. 
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Fic. 4. Spanwise distribution of section lift coefficient. 


Fig. 3. A.R. = 3.45, 46.4° sweep, taper ratio 2.4 to 1. Excel- 
lent agreement of the method of Eq. (1) with Weissinger and 
experiment is noted over the inboard 70 per cent span. Over the 
outboard portion, the method of Eq. (1) is slightly closer to 
experiment than Weissinger. 

Fig. 4. A.R. = 10, 35° sweep, taper ratio 3:1. The method 
of Eq. (1) compares well except at the root, where the method 
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U = velocity of piston, ft. per sec. 
bee 425° 96 3° —-- SHRENK (EXTENDED) Po = initial static pressure, Ibs. per sq. ft. 
Pipers — —EXPERIMENT( @ 16 $*,C, +0 585) po = density, slugs per cu.ft. 
a" /N 8 = angle between the tangent to the surface of the airfoil and the chord 
« ? a a" f- ~\ of the airfoil 
Lae ‘\ M = Mach angle 
40 a 7 : , P ro 
, Cae \ q = dynamic pressure = (1/2)p,V? 
oe _ x v = volume, cuft. 
\ 
06 ' Fig. 1 illustrates schematically a thin wedge-type airfoil in two 
nh \ dimensions moving at supersonic velocity (V). In position ‘‘A” 
\ es (eae 
os \ two inclined waves, the Mach waves, are observed originating 
1 . , 7 4 1. 1 ‘ at the airfoil nose. The Mach waves are envelopes of the dis- 
‘ “4 = ss = =e id Tr a '°  turbances produced by the nose of the airfoil as it travels for- 
Fic. 5. Spanwise distribution of section lift coefficient ward; these waves make an angle M, the Mach angle, with the 


of Theilheimer rises asymptotically and must be faired. No 
experimental data are available. 

Fig. 5. Delta wing, A.R. = 
with experiment at all points out to 90 per cent span at low values 
of C,. Agreement is actually better than the figure might indi- 
cate, because the Schrenk (extended) curve should be raised 
slightly to give the same value of Cy, as the experimental curve. 
The Schrenk curve is faired to zero from the 90 per cent span sta- 


2, 56.3° sweep. Good agreement 


tion. 

The authors wish to note that sufficient time has not been 
available to make complete calculations at many aspect ratios, 
taper ratios, and twists; hence, Eq. (1) must be used with caution 
The authors would welcome comments. 
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An Approach to Supersonic Airfoil Theory 


W. E. Strohmeyer and D. R. Gero 
North American Aviation, Inc., Los Angeles, Calif., and Goodyear 
Aircraft Corporation, Akron, Ohio, respectively 


April 13, 1949 


HIS NOTE PRESENTS an approach to supersonic airfoil theory. 

The basic equation, as commonly associated with Ackeret, 
used in determining the aerodynamic forces on thin airfoils 
traveling at supersonic velocities, is obtained by making use of 
the fundamental energy equation. 

The purpose is to develop the basic expression for the aero- 
dynamic forces of thin airfoils traveling at supersonic speeds 
utilizing a concept that provides a convenient picture of the 
mechanism by which these forces are developed. The results are 
similar to those developed by Taylor! and commonly associated 
with Ackeret. It is believed that this approach not only results 
in a derivation that is more easily understood but also provides 
a greater appreciation of the assumptions and limitations of the 


theory. 
SYMBOLS 
V = supersonic velocity, ft. per sec. 
T = time, sec 


local velocity of sound, ft. per sec. 


direction of motion of the airfoil. This angle may readily be 
determined by noting that at position 
propagated radially from the airfoil nose at a velocity equal to 
the local velocity of sound. At an increment of time later AT, 
during which the airfoil has traveled a distance (ATV) to posi- 
tion ‘‘B,’’ the wave front has progressed to point A’ covering a 
distance (ATC). 
nose disturbances, line BA’ is the Mach wave, and the Mach 
angle is therefore the angle between line BA’ and the direction of 
flight AB. It will also be noted that, at position B, the surface 
of the airfoil has succeeded in displacing the fluid from point 
“A” toa’. This fact is the key to further analysis of the problem 
since the pressure arising from this displacement acts on the sur- 
face to produce the aerodynamic forces. 

In the preceding paragraph, it was indicated that the pressure 
acting on the surface of the airfoil is caused by the air being dis- 
placed in the general direction of the advancing sound wave. To 
examine this mechanism in detail, it is convenient to use an anal- 
ogy whereby the surface of the airfoil is replaced by a piston act- 
ing upon a column of air of unit cross section as shown in Fig. 2. 

If the piston is made to follow a velocity step function so that 
the velocity varies instantaneously from zero to a small value U, 
a pressure wave will be propagated forward at a velocity equal 


‘‘A”’ pressure waves are 


Since the Mach wave is the envelope of the 
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FIG 4 
PISTON ANALOGY 


to the local velocity of sound (C). If, after this initial disturbance, 
the piston is moved forward at the constant velocity U, the pres- 
sure of the air enclosed between the piston face and the pressure 
wave front will be increased from Py to Py + AP. This can be 
seen by noting that, in time AJ, the pressure wave has advanced 
a distance (ATC) affecting a like volume of air—i.e. (distance 
traveled times area) = ATC, and during this same interval the 
piston has advanced a distance (UAT) compressing the original 
volume (ATC) by an amount (UAT). Since the “affected”’ vol- 
ume and the portion compressed are both linear functions of time, 
the per cent compression remains constant if velocities C and U 
are constant. If no further compression of the gas occurs behind 
the sound wave, it is evident that the compressed gas must be 
moving with the same velocity as the piston. 

It is now possible to write an energy expression equating the 
work done by the piston to the gain in the kinetic plus the 
increase in potential or stored energy of the gas. 

piston work = A kinetic energy + A potential energy 
The work done by the piston per unit time is 
Force X velocity = U(P» + AP) 
In unit time, a mass of air equivalent to that affected by the 
pressure wave (poC) is given a velocity equal to that of the piston. 
Therefore, the gain in kinetic energy per unit time equals 
A K.E. = poCU?/2 

The stored energy or the increase in potential energy is that 
due to the compression of the gas. If the degree of compression 
is small, the pressure may be considered a linear function of the 
volume, as shown in Fig. 3—i.e., 

AP = kAo 

The increase in stored erergy or potential energy is equivalent 

to the work required to compress the gas; for a unit area piston, 
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it is the product of the average pressure on the piston and the 
distance traveled. 


AP ; 
AP.E. = (Po+ =) AV, AV = Ude 


A P.E. = PoU dt + (APU dt/2) 


Hence, the increase in stored energy per unit time = Pol’ + 
(A PU/2). The complete energy equation then becomes 
, ' 0CU? APU 
PoU + APU = ™—— + Pl’ + 
or 
AP = poCU (1) 


Thus, the pressure acting on a surface behind a pressure wave 
depends upon the initial density po, the velocity of wave propaga- 
tion (C), and the velocity of the surface normal to the advancing 
pressure wave (U). It must be remembered that this solution 
assumes that the pressure is a linear function of the volume which 
is valid only for small changes and thus limits the velocity (U) 
to small values compared to velocity (C). 

Fig. 4 represents a wedge airfoil traveling at supersonic velocity 
V at zero degrees angle of attack with the accompanying Mach 
wave. The airfoil is shown in two positions illustrating the man- 
ner in which the surface acts as the piston described above. Lines 
“‘aa’’ represent the outline of a hypothetical cylinder extending 
to infinity, the Mach wave represents the pressure wave traveling 
at a velocity (C), and the surface of the airfoil is analogous to the 
piston. The velocity (U) of the surface is found from Fig. 4 to be 


U = V sin B/cos (M — 8) (2) 
For thin airfoils sin 8 = 8, and, if the Mach angle is extremely 
large compared to 8, it is permissible to write, 


U = V£B/cos M (2a) 


where 1/cos M = V/V v2 =-@ (see Fig. 4). 
Substituting for (u) in Eq. (1), 


AP = 298/*/(V2/C) — 1 (3) 
which is the Ackeret formula.? 
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Errata—Helicopter Control and Stability in 
Hovering Flight 


R. H. Miller 
Associate Professor of Aeronautical Engineering, Massachusetts 
Institute of Technology, Cambridge, Mass. 


May 12, 1949 


Au OF ERRORS in the article ‘“Helicopter Control and Stabil- 
ity in Hovering Flight,”” by R. H. Miller (JOURNAL OF THE 
AERONAUTICAL SCIENCES, Vol. 15, No. 8, p. 453, August, 1948) 
follows: 

Page 459: (1) Equation for Cg: ox should read Ar. (2) J, 
and J, do not include the moment of inertia of the blades but must 
include the weight of the blades, assumed to act at the rotor hub, 

Page 460: (1) In the definition of w, after Eq. (7), Q should 
read 2. (2) Eq. (8), u? should read: U?. 

Page 461: (1) Eq. (12), the second to the last term, — CQ, 
should read —Cdpr. (2) Eqs. (13a) and (13b), the second terms 
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should be divided and not multiplied by 2. (3) Eqs. (14) and 
(15), the denominator of the last term should read 2(1 — £)%p 
and not 2(1 — &)lp. 

Page 462: (1) Integral for P, Insert R before dr. (2) In the 
expressions for Ho, Hi, etc., A should read Ar. (3) In the expres- 
sion for Hy, Cp,(1 — &) should read (Cp,/a)(1 — &). (4) Eq. 


Al AI\? 
(16), should read 2 
I 7 


Page 463: (1) In the expressions for Ax: multiply the first term, 
v2, by [1 + m(e/d)] and replace Ms by [Ms + 2m(e/d)]. With 
this modification the qualification ‘‘neglecting the last two terms 
of Eqs. (14) and (15)... .’’ may be eliminated. (2) In the ex- 
pression for Rx, ¢ should read ¢ + w,? + 6B. 

Page 464: (1) In the expression for J, M4 should read Mg. 

Page 465: (1) In the stability determinant, 7rCp,/4 should 
read rCp,/4a. 

Page 466: 
read £(a;). 


(1) In text, after “‘Laplace transform,’ £(a;) should 
(2) The expressions for G(v) and g(v) should read 
G(v) = Xv + (HX — JZ) 


g(v) = Yr + (HY — JV) 


where 


ae, H, D.C), He 
~ YS 14 /A)LA AB 1 + (¢/A)? 


sD i BoC } 
A* A 


Page 467: (1) The statement after the simplified expressions 
for XY and Y should read: ‘In determining initial response, 
H and J may always be assumed zero.”’ (2) Eq. (20), X and Y 
should be multiplied by 2?. 

Page 468: (1) Equation before Eq. (22), the fourth \ should be 
squared. (2) Eq. (22), —[Q°kX — JH,Q'] should read —[Q?kX 
— JH]. (3) In the expression for D’, multiply 2(@) — 
5Bo) by rR. 

Page 469: 
out. (2) Following Eq. 


(1) Eq. (24), f should be multiplied by 2 through 
(24), the first roots given should be 
\ = — 0.585 + 3.537, ete. (3) First column, fifth line from the 
bottom, a;/@; = 4.5, not 1.5. (4) In the expressions for (a; — 81) 
and fp, the terms containing C,,, should not contain a. 

Page 470: 
should read Eq. (25). 

Page 471: 
Cy. (2) In expression for a/6,, u should read x. 
lowing Eq. (2A) “‘right’’ should read “‘left.” 
(1) pb/2u should read pb/2U. 


(1) First column, fourth line from bottom, Eq. (23) 


(1) In expressions for m, and zg, C; should read 
(3) In line fol- 


Page 472: (2) ¢/6, is negative 


Concerning Linearized Supersonic Flow Solu- 
tions for Rotationally Symmetric Bodies 


W. H. Dorrance 
Research Asscciate in Aerodynamics, Aeronautical Research 
» Center, University of Michigan, Ann Arbor, Mich. 


May 19, 1949 


UMEROUS AUTHORS!~* have formulated solutions to the lin- 
N earized partial differential equation of compressible poten- 
tial flow applicable to rotationally symmetrical bodies when the 
incident Mach Number exceeds 1. The solution of von Karman 
and Moore for such bodies aligned at zero angle of attack, and 
the Tsien? or Ferrari’ solution for bodies aligned at a small angle 
of attack with the incident flow can be shown to have a common 
origin in certain solutions appearing in the theory of wave 


mechanics. Unfortunately, the Tsien or Ferrari solution has 
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been misrepresented to some extent in certain widely distributed 
texts.5,® The purpose of this note is to clarify this misrepresen- 
tation by examining the origin of the solutions formulated by the 
above-mentioned authors. 

The linearized partial differential potential equation governing 
these solutions in cylindrical coordinates is 


m 0° 1 0 


o% 1 oF 
(M? — 1) ht = ; . 
Ox? or? ror 


(1) 
r? 06? ‘ 


where M = free stream Mach Number and (x, 7, @) are cylindrical 
coordinates. 

For the axially symmetric flow about a rotationally symmetric 
body, the flow is independent of @ and Eq. (1) becomes 


) O*h0 a O*h0 1 Odo 
ng r or 


(M?—-1 = ; 
Ox? or? 


von Karman and Moore employed the simple substitution ¢ = 
x/V to change Eq. (2) into the familiar form of the two-dimen- 
sional wave equation for infinitesimal disturbances. They formu- 


lated the solution below arising from a solution to this wave 


l 0 
f f(x — Br cosh u)du (3) 
4dr cosh~! (x/Br) 


where 6 = VM? — 1. 

If the substitution = x — £r cosh u is used in Eq. (3) this 
solution takes the form of the source-sink solution, appearing in 
incompressible flow solutions. 


x — Br . 
iad mf SO 
0 V (e — &)* — 6%? 


This solution represents the potential at a field point (x, 7) due 
to a distribution of sources and sinks along the x-axis at points 
£ between 0 and x — @r. 

Following in form the procedure outlined by Lamb,’ the general 
solution to Eq. (1) when the Mach Number exceeds 1 is 


equation. 


go = 


n 


z. d.r* cos sO + yr* sin sO ( 


s=1 


or 


= 


where ¢, and y, must satisfy the equation below obtained by 
substituting Eq. (5) into Eq. (1). 


o? ae L  2s+12 4 t . 
ate vos = Pats Ve + © ore a | (6 


A known solution to Eq. (6) is 


ds = [(1/r)/(0/Or)]"do 


(7) 
where ¢po is identically the solution to Eq. (2) given by von Ka4r- 
m&n and Moore. 
By taking the index s equal to 1 in Eq. (5), Eq. (7) yields the 
solution below formulated by Tsien and Ferrari 
@ = (0¢0/07r) cos 6 (8) 


The differentiation of ¢o is performed on the integral in the form 
of Eq. (3) because the integral in Eq. (4) is improper since the 
denominator vanishes at the upper limit. The solution to Eq 


(1) is then 


B cos 6 0 ” 
—— f'(x — Brcosh u) cosh u du (9) 
4 cosh~! (x/fr) 


if f (0) can be taken as O, a reasonable assumption since these 
solutions within the rigor of the linearization are restricted to 


o¢= 


sharp pointed bodies. 
Calling —& = x — Br cosh u and f’(€) = g(é), solution (9) be- 


comes 


x — Br 
‘iis —— 
@ = (cos 0/4rr) (t)(x — &)dt/V (x — £)? — 8? 
vO 
10 


ce 


a 


tributed 
epresen- 
d by the 


verning 


(1) 
lindrical 


mmetric 


bo 


ion? = 
-dimen- 
’ formu- 
is wave 


(3) this 
iring in 
(4) 


,r) due 
points 


general 


1 K4r- 


ds the 


(8) 


» form 
ce the 
Oo Eq. 


these 
ed to 


9) be- 


10 
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This solution has been incorrectly represented in references 5 
and 6 as 
x — Br 


g(t)dé/[(x — &)? — pr]*/* (11) 


¢ = (6°r cos 6/4) 
0 
The integrand of this incorrect solution resembles in form the 
integrand of the subsonic doublet solution obtained by differ- 
entiating the source-sink solution integral (4) with respect to r 
under constant limits. When the Mach Number is less than 1, 
this procedure is valid because the limits of integration are con- 
stant for a body of finite extent. When Mach Number exceeds 
1, the upper limit becomes x — §8r as a consequence of Eq. (2) 
changing from the elliptic type to the hyperbolic type characteris- 
tic of supersonic flow fields. Only when Mach Number exceeds 
1 can solutions of the two-dimensional wave equation be adapted 
to Eqs. (1) and (2) in the manner discussed herein. As such, 
solution (10) is the proper form of the solution to Eq. (1) when 
Mach Number exceeds 1 and should be used instead of the form 
of Eq. (11) whenever Eq. (11) appears. 
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Quasi-Stationary Airfoil Theory in 
Compressible Flow 


John W. Miles 


Department of Engineering, University of California at Los Angeles 
April 29, 1949 


[ WAS RECENTLY POINTED OUT that the use of steady flow the- 
ory could lead to erroneous results in the calculation of sta- 
bility derivatives, and the correct results for quasi-stationary in- 
compressible flow over a thin airfoil were obtained as limiting 
results of the previously known flutter results.1 These results 
have now been extended to the case of subsonic, compressible 
flow by obtaining a solution to the Possio integral equation,? 
considering only terms which are first order in frequency. 

The lift and quarter-chord moment coefficients on a thin air- 
foil with leading edge at x = —1 and trailing edge at x = +1in 
a subsonic flow of velocity U and Mach Number / with down- 
Wash distribution w(x) exp (twt) are given by 


CG = 2(1 — M2)—"/s Sa — cos ¢) + 
0 \ 


ik(1 — M?)-1 sin? ¢ + ik(1 — M?)-'(1 — cos oly + * oh 


(— cos ¢) 
ios) “ ran | t ae ¥' de + O(k* log &) (1) 





. cludes a term that is zero order in frequency. 
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1 7 
C= re — M%)-'/: [(cos ¢ — cos 2y) — ik(1 — M*)-!X 


0 
— cos ¢) 
(1 — cos ¢) sin? ¢] = = = dg + O(k? log k) (2) 
2(1 — M?) 
F(M) = M? + toe | — (1 — M?)? log, X 
[? + (1 - | 
a (3) 
M 
k = we/2U (4) 


where & is the reduced frequency, based on the semichord (c/2), 
and y is Euler’s constant. These results reduce to the well- 
known Munk formulas,’ together with the Prandtl-Glauert 
factor (1 — M?)—'/3, for k = 0. 

For the important case of pitching about the quarter chord— 


W(x) = Ul1 + [(1/2) + (2x/c)] ik} (5) 
Eqs. (1) and (2) reduce to 


Ca = 2r(1 — a)-") (1 + ik) + ik(1 — M%)-! x 


1 k oe 
= + vy + logd =} — F(M) f (6) 


= —(irk/4)(1 — M2)—*/(2 — M2) (7) 


Cre 


For plunging of the airfoil—i.e., 


w,(x) = ikU (8) 

the results are 
Cy = 2eik(1 — M?)—"/2 (9) 
Cm, = 0 (10) 


In interpreting these results, it may be remarked that nonsta- 
tionary effects need be considered only when the downwash in- 
In this case, how- 
ever, the term introduced is logarithmic and will be increasingly 
important for small k, corresponding to high-speed flight. More- 
over, the compressibility correction of this term is approximately 
a- M?)—*/2 rather than (1 — M?)~—'/2. It follows that the ef- 
fects under consideration will be particularly important at high 
speeds. 

As an example, consider the calculation of the damping of a 
rotary motion of a tail surface about a center five chord lengths 
ahead of its quarter-chord. Carrying out the calculations with 
the aid of the above results and also with the results obtained by 
setting k = 0 in Eqs. (1) and (2) and designating the damping 
derivatives as C,,, and C°,,,, respectively, it is found that the latter 
calculation overestitnates the damping by a considerable margin. 
Typical numbers are 





M* | 0.1 | 0.01 
C, a ae | > = ; | > } 
- = 0 } 0.828 | 0.618 | (11) 
Con ee ts ct oe 
° a 
|} 0.7 | 0.342 


0.662 


While these ratios neglect induction effects and wing interference, 
they are indicative of the errors that may be expected when 
steady flow theory is used for the calculation of tail stability 
derivatives. Similar results may be expected for the damping 
in pitch of a swept wing. 
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First Order Wing-Body Interference Effects 


E. V. Laitone 

Associate Professor, Mechanical Engineering, University of 
- California at Berkeley 

April, 1949 


a py PROCEDURE GIVEN BY Ferrari! on wing-body interference 
at supersonic speeds is not valid for rectangular wings pro- 
jecting out less than a body diameter (a common configuration 
for a supersonic missile). This is so because when the wing span 
approaches the body Ferrari's assumption ‘‘that the induced 
field generated by the wing is that which would exist around the 
wing if it were placed in the uniform stream alone’’ cannot be 
true, so the interference of the wing on the body, when computed 
on this basis, will be too large. Also, whenever the Mach cones 
from the wing tip intersect downstream on the body, the resulting 
wing downwash acting on the body will also decrease the inter- 
ference corrections. Consequently, for the case where the total 
wing span is only slightly larger than the body diameter, it is 
probably best to use the simple first order correction as obtained 
from the additional mass? and linearized theory.* 
Then with the notation in Fig U’ is the flow in- 
duced by the wing we have: 

(1) Influence of the Normal Velocity Induced by the Wing on 
the Body.—For either subsonic or supersonic flow the first order 
expression for the transverse force distribution is the same as for 
incompressible flow (see Eqs. 44 and 58 of reference 3). There 
fore the apparent mass concept gives the interference, based on 


the wing area S and the wing chord C, from reference 2 for any 


apparent 
1, where e = w 


wing-body combination 


AiC, = —2e(L)[A(L)/S] 


: A(L)(L G 2 L 
AiCug = 2e(L) — _— — i — e(x)A(x)dx 
Ss \¢ ¢ SC 0 
where A = ar? is the cross-sectional area of the body. Now for the 
special case of supersonic flow with a rectangular wing mounted 
on the cylindrical afterbody as in Fig. 1, we have 
L 


e(x)A(x)dx = aCA(L) 


0 
provided the wing tip Mach cones do not intersect on the body so 
e = O everywhere ahead or behind the wing root Mach cones 
within which e = a. 

If the body extends a long distance behind the rectangular 
wing, then in supersonic flow with a uniform lift distribution 
everywhere (except in the small triangular regions formed by the 
wing tip Mach cones which may be neglected for Mach numbers 
sufficiently greater than unity) 
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e(L) + e(@) = C,/xA.R. 


However, for subsonic flow the lift distribution is more near}, 
elliptic, so for a long body 


e(L) = e(~) = 2C_,/rA.R. 


(II) Influence of the Longitudinal Velocities Induced by the 
Wing on the Body.—For an infinitely long span rectangular wing 
on a body of small diameter, the lift increment is given by th. 
two-dimensional lift distribution acting on the actual area coy 
ered by the body; so for supersonic flow with b and C > D 


AoC, — (da B)(DC, S) 


provided the body extends sufficiently behind the wing so as to 
include all the Mach cones emanating from the wing-body junc- 
ture. However, for a finite wing span, A,Cz must be modified 
so for the limiting case of a small span rectangular wing we can 


write 
: ta DC b-—D 
Aol L = * — — 
BS ~/(b — D)? + D? 
AcCna = — (Gy C) AsCr 


where G,, is the distance from the center of pressure of the wing 
(only a small distance ahead of the mid-chord position) to the 
center of gravity (G). 

(III) Influence of the Normal Velocities Induced by the Body 
on the Wing.—As shown in reference 2, the normal velocities in- 
duced by a long slender body may be taken as those correspond- 
ing to two-dimensional flow about a circular cross section, so that 
for subsonic or supersonic flow on any wing plan form 

b/2 
C.(D/2y)*¢ dy 
D/2 


A3CL = 2/S 


Then for supersonic flow on a rectangular wing, with Mach Num- 
bers sufficiently greater than unity so the tip Mach cones affect 
only a negligible region, 


ta Di(b— D)C 


A3sC, = 
_ Bb Bs 


AsCmg =~ — (Gy “oy A;CL 


(IV) Influence of the Longitudinal Velocities Induced by the 
Body on the Wing.—As shown in reference 3 the longitudinal 
velocity of the body is of the order r? log 7, so it may be neglected 
in the linearized wing theory so that 


AyCr =Q= AiCng 


The total lift and moment based on these corrections must then 
be computed by considering only the external wing area actually 
projecting out from the body. For example, the total lift would 
be given for the configuration shown in Fig. 1 as 


ae da 1 4 a (6 — D)C - A(L) + AC. + 
= 2BA.R. S ~~ s ae 
AsC, + AsC 


It is important to note that the aspect ratio should be taken 
as 


A.R. = (b — D)/C 


because as b — D the body acts as an infinite plane of reflection 
while for 6 > D the wing acts independently of the body 

The above first order interference corrections are also applicable 
to large span wings provided they remain entirely within the 
Mach cones emanating from the nose of the body. If the wing 
extends beyond the body’s Mach cones, then the upper limit of 
integration for A;C;, must be correspondingly reduced. 
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On the Oscillating Aileron at Supersonic 
Speeds 


John W. Miles* 
Department of Engineering, University of California at Los Angeles 


May 3, 1949 


HE SOLUTION OF THE oscillating, supersonic, rectangular 
ae problem! was recently announced. As a continuation 
of this investigation, a solution has been carried out for the os- 
cillating rectangular control surface that has its outboard edge 
adjacent to the free stream and its inboard edge joined to a wing 
through a sealed gap. 

The calculation of the pressure distribution in the zone of ac- 
tion of the outboard edge is identical with the wing edge calcula- 
tion, while the calculation in the zone of action of the inboard 
edge is a problem of the first kind, in the sense that the downwash 
within the zone of influence of every point within the latter zone 
of action is prescribed. In calculating the lift and moments on 
the wing due to the deflection of the control surface, we must in- 
clude the pressure acting on the wing in the zone of action of the 
inboard edge, but this pressure does not contribute to the hinge 
moment, 

The results may be stated as follows: 

(1) The force acting on a spanwise strip bounded by the two 
Mach lines from the leading, inboard corner is just equal to the 
force that two-dimensional (‘‘strip’’) theory would predict for 
that half of this strip lying on the control surface proper. 

(2) The force acting on a spanwise strip bounded by the in- 
board edge and the inboard Mach line from the leading-edge cor- 
ner requires an aspect ratio correction which is (2/m) times the 
correction that would be applied if the edge were outboard, the 
latter correction being given in reference 1. 

(3) The lift coefficient, referred to the product of the dynamic 
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pressure and the control surface area, but including the lift in 
the entire zone of action of the control surface, due to a unit 
(amplitude) rotation of the control surface about its leading edge, 
is given by 
Cig = 4(M* — 1)~* ff [go(2xu) — (4xA.Rg)-? X 
gi(2xu)) filk, u)du 
(4) The coefficient for the moment about the leading edge 
defined as in (3) with the control surface chord as the reference 
moment arm, is given by 
ont ‘ = 
Cus = —4(M? — 1) hd ,% [go(2xu) = (4xA.R.,) ae 
gi(2xu)] fo(k, u)du 
(5) The coefficient for the rolling moment about the control 
surface mid-chord line, defined as in (3) with the control surface 
span as the reference moment arm, is given by 
Cig = (2xA.R..)-"(M? — 1) =a So’ (gi(2xu) — 
(4xA.R.,)~ !go(2xu)] fi(k, u)du 
(6) The hinge moment coefficient, defined as in (4) but includ- 
ing only the forces on the control surface proper, is given by 


Cop = —4(M? — 1)— 7? J)" [go(2xu) — [1 + (2/x)] X 
(4xA.R.,)~! gi (2nu)] filk, u)du 

k = (w/u) X (semichord) 

K = M(M? — 1)—%® 

A.R. = (M?— 1) AR. 

Silk, u) = (1 + 4tk — 2k?) + 4(Rk? — ik)u — 2k2u? 

folk, u) = 2[tk — (2k?/3)] + (1 + 2k*)u — 2iku? — (2/3)k*u3 

go(.X) = exp (—iMx) Jo (x) 

g(x) = exp (—7Mx) sinx 

g(x) = x exp (—iMx) Ji (x) 


where M is the Mach Number, & is the reduced frequency, w is 
the angular frequency (implying the time dependence exp(iw/), 
A.R. is the true aspect ratio of the control surface, and J,,(x) is 
Bessel’s function of order n. 

The integrals involving go arise in the two-dimensional problem 
of the oscillating airfoil and are therefore known.? The integrals 
involving g2 may be reduced to integrals of go after integration by 
parts, and the integrals involving g; are elementary. 
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THE 


JOURNAL OF 
SUGGESTIONS for CONTRIBUTORS 
to the Publications of the 
INSTITUTE of the AERONAUTICAL SCIENCES 


AERONAUTICAL 


SCIENCES—AUGUST, 1949 


The Institute of the Aeronautical Sciences invites both members and nonmembers from any 
country to submit papers for publication in the JOURNAL OF THE AERONAUTICAL SCIENCES and the 
AERONAUTICAL ENGINEERING REVIEW. The Institute, following the practice of other societies, 


does not pay for contributions. 


The following directions for the preparation of papers, if followed by authors, will save 
correspondence, avoid the return of papers for changes, minimize the work of preparation for the 
printer, and save the expense due to the charges made for ‘‘author’s corrections.”’ 





Manuscripts: The original typewritten copy of the paper 
desired, double or triple spaced on one side of white paper sheets, 
consecutively numbered. There should be wide margins to 
allow for the marking of directions to the printer. Correcting, 
changing, or adding to papers after they are in type is costly. It 
is, therefore, imperative that papers submitted be in final form. 
Typographical errors may be corrected on proofs, but if authors 
wish to add material, they may do so at their own expense. In 
mailing, drawings may be rolled, but manuscripts should be sent 
flat. Send by first class mail (register, if you wish, for your own 
protection) to the Editorial Office, Institute of the Aeronautical 
Sciences, 2 East 64th Street, New York 21, N.Y. All manu- 
scripts will be examined by the Editorial Committee and by the 
Editor. Authors will be advised as promptly as possible whether 
the paper is acceptable for publication. 


TitLes: The title of the paper should be brief. The name 
and initials of the author should be written as he prefers. The use 
of the full name of an author is advocated because of the fre- 
quent duplication of initials and surnames which sometimes 
makes it difficult to establish the identity of the author. This is 
particularly important for large annual indexing and abstracting 
services. All titles and degrees or honors are omitted. The 
name of the organization with which the author is associated 
should be placed after his name on a separate line. The date on 
which the paper is received will be inserted by the Editor. The 
author’s title should be indicated in a footnote. 


SUMMARIES OR AsBsTRACTS: An abstract to be printed at the 
beginning should accompany each article. It should be ade- 
quate as an index and asa summary. It should contain a state- 
ment of major conclusions reached, since summaries in many 
cases constitute the only source of information used in com- 
piling scientific reference indexes. Abstracts printed in other 
journals, especially foreign, in most cases, consist of summaries 
from printed papers. The summary should explain as adequately 
as possible the major conclusions to a nonspecialist in the subject 
and should contain from 100 to 300 words, depending on the 
length of the paper. 


SUBHEADINGS: Subheadings should be inserted by the author 
at frequent intervals. The work of editorial preparation will be 
simplified by the author providing many subheadings. 


MATTER USUALLY DELETED: Photographs or illustrations of 
little technical interest and not showing advances in general prac- 
tice. Too detailed tabular matter (general results of such tables 
may be included in the text). Lengthy descriptions of materials 
or processes or of preliminary experiments or theories that pre- 
ceded final results; salient features only are of interest. 


REFERENCES AND FOOTNOTES: References should appear as 
footnotes only, numbered consecutively, grouped together at the 
end of the manuscript. The arrangement should be as follows: 
(for books)—! Durand, W. F., Aerodynamic Theory, 1st Ed., 
Vol. 1, p. 23; Julius Springer, Berlin, 1934. (For magazines)— 
1 Englund, C. R., Crawford, A. B., and Mumford, W. W. Some 
Results of a Study of Ultra-Short-Wave Transmission Phenomenon, 
Proc. I. R. E., Vol. 20, No. 12, pp. 481 and 482, March, 1933. 
Please give author, title, edition, volume, page, publisher, and 


date of publication as indicated. Omission of one required fact 4 
causes much extra editorial work and possible inaccuracies. 


ILLUSTRATIONS: Illustrations should accompany manuscripts 
and each should always be referred to in the text, by number. 
Drawings or graphs should not be larger than 12 X 16 inches, and 7 
must be made with jet black India Ink on white paper or tracing 
cloth, the latter being preferred. Do not use typewriter for 
lettering. The smallest lettering on 8 X 10 inch figures should 
be no less than '/, inch high. Cross-section paper (white with 
black lines) may be used, but should not have more than 4 lines 
per inch. If finer ruled paper is used, the major division lines 
should be drawn in with black ink, omitting the finer divisions, 
In the case of finely ruled paper, only blue-lined paper can be 
accepted. Tracing paper and blueprints are not acceptable. 
Lettering and all markings must be large enough to be readable 
after reduction. Méail rolled or flat, never fold. Drawings that 
cannot be reproduced (including pencil drawings) will be returned 
to the author for redrawing, thus delaying publication of the 
paper. Photographs should be distinct and show clear black and 
white contrasts. They must be on glossy white paper. Avoid 
round and oval photographs. 


CAPTIONS AND LEGENDS: Legends or captions must accom- 
pany each drawing or photograph submitted. If written on the 
drawing or photograph, they should be placed below and well out- 
side the part to be reproduced. Each table should have a caption 
such as Table 1, Table 2; Table 3, etc. Captions should be com- 
plete in themselves so as to make the data intelligible to the 
reader without reference to the text. A duplicate list of captions 
for figures should be included as the last page of the manuscript. 
Use ‘‘Fig. 1’’ (not Figure 1), Figs. 3 and 4, etc., in both the text 
and the numbering of illustrations. In the text, ‘‘Eq. (1),” 
or “Eqs. (1) and (2)”’ are preferable to ‘‘Equation (1).’’ In cap- 
tions and legends, except for ‘‘Fig.”’ and ‘‘Eq.,’’ and in table head- 
ings, write all words in full; do not abbreviate. 


MATHEMATICAL WorK: Only the simplest formulas should 
be typewritten, all others should be carefully written in pen 
and ink, the writing to be large enough so that ample room is 
provided to mark mathematical matter for the printer. A con- 
siderable space for marking should be allowed above and below 
all equations. All complicated equations should be repeated on ~ 
separate sheets with plenty of space left for marking. The solidus 
should be used for simple fractions appearing within the text. 
Make all expressions clear to the typesetter. Greek letters used 
in formulas should be clearly designated by name on the margin 
of the manuscript. All symbols should be clearly written and 
carefully checked. The difference between capital and lower- 
case letters should be clearly distinguished and care taken to 7 
avoid confusion between zero (0) and the letter (0), between the ~ 
numeral (one) and the letter (ell) and the prime (’), between alpha 
and a, kappa and k, u and mu, v and nu, n and eta. All sub- 
scripts and exponents should be clearly marked and dots and bars 
over letters or mathematical expressions should be avoided. 
Avoid complicated exponents and subscripts. When it is neces- 
sary to repeat a complicated expression, it should be represent 
by some convenient symbol. F 


NOMENCLATURE AND ABBREVIATIONS: Standard abbreviations 
should be used, and it should be noted that most abbreviations 
are lower case such as m.p.h., b.m.e.p., i-hp., b.hp., hp., . . . ete. 





ed fact 


iscripts 
lumber 
1es, and 
tracing 
iter for 
should 
te with 
4 lines 
yn lines 
visions, 
can be 
ptable 
adable 
gs that 
turned 


of the 
ck and 
Avoid 


1ccom- 
on the 
ell out- 
aption 
e com- 
to the 
iptions 
script. 


should 


yom is 
A con- 
below 
ted on” 
olidus 
> text. 
s used 
nargin 
n and 
lower- 
cen to 
on the 
alpha 
1 sub- 
1 bars 
oided. 
neces- 
ented 


itions 
tions 
etc. 





